I am solving system of two coupled equations from fluid dynamics.
$p_0 p_0'=-\dfrac{32 \beta}{R^4}$
$(p_0 p_1)'=-\dfrac{2-\sigma_v}{\sigma_v}p_0'\dfrac{8}{R}$
where $\beta$, $\sigma_v$ are constants. $R=R(z)$ radius dependent on $z$-it is not constant, where $z$ is longitudinal coordinate and $R=r-z(r-1)$, where $r$ is inlet radius. All these equations are in dimensionless form.
As a solution I need to have pressures $p_0$ and $p_1$, these pressures are first and second approximations of pressure, I got it with series expansion of pressure. And I have boundary conditions, that at the inlet $p_{0}|_{z=0}=p_i$ and at the outlet $p_0|_z=1$, for second approximation $p_1|_{z=0}=p_{1}|_{z=1}=0$. Because of this definition I implemented shooting method and got neccessary conditions.
I tried to use MATLAB code from address https://www.mathworks.com/help/symbolic/solve-differential-algebraic-equations.html just instead of ode15 which is used on link, I need to implement ode45 (because I need to implement Runge Kutta method).
**
My problem
** is because ode45 is everywhere in matlab implemented with variables which depend on time $t$, and my pressures $p_0$ and $p_1$ are only functions on longitudinal coordinate $z$. Do I need in matlab code write my variables as $p_0(z)$ or just $p_0$, and do I need to change varible $z$ in my equations with time $t$ because matlab ode45 works only with variable $t$? Does diff in matlab accept that derivation is not according to time $t$?
I tried couple variants, and every time I got error.
Do you have any other helpful advice for solving this system of ODE with Runge Kutta method?
My code is here, I believe there is something wrong with ode45 part of the code, but I cannot conclude what.
> syms R(z) z r beta sig_v p0(z) p1(z);
>
> R(z)=r-z*(r-1);
>
> eqn1=p0(z)*diff(p0(z),z,1)==-32*beta/(R(z)^4);
> eqn2=diff(p0(z)*p1(z),z,1)==-(8*(2-sig_v)*diff(p0(z),z,1))/(sig_v*(R(z)));
> eqns=[eqn1 eqn2];
>
> vars = [p0(z); p1(z)]; origVars = length(vars);
>
> isLowIndexDAE(eqns,vars) [DAEs,DAEvars] = reduceDAEIndex(eqns,vars)
> isLowIndexDAE(DAEs,DAEvars) pDAEs = symvar(DAEs);
> pDAEvars =symvar(DAEvars);
> extraParams = setdiff(pDAEs, pDAEvars)
> %index is low, so we can go further
>
> f = daeFunction(DAEs, DAEvars, beta, r, sig_v);
> beta=1; sig_v=1; r=2;
>
> %%%%%%%%%%%% boudary condition for p and p'
> solinit = bvpinit([0 1],[0 1]);
> sol = bvp4c(@twoode,@twobc,solinit);
> xint = linspace(0,1);
> yint = deval(sol,xint);
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> % ode 45 solving
> options=odeset('RelTol', 1e-8, 'AbsTol', [1e-8 1e-8]);
> g=@(z,p0(z),p1(z), beta, R, sig_v)[eqn1; eqn2];
> [p0(z),p1(z),z]=ode45(@(z,p0(z),p1(z)) g(z,p0(z),p1(z), beta, R, sig_v), [0 1.5],[1 1/2]);