Two couple equation first order Runge Kutta Matlab troubles?

Question

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]);

Answer 1

You can integrate the first equation directly, $$ (p_0(z)^2)'=-\frac{64}{(2-z)^4}\implies p_0(z)^2=p_0(0)^2+\frac{64}{3(2-z)^3}-\frac83 $$ As you can see, you only have one integration constant free that you can use to either fix $p_0(0)$ or $p_0(1)$, but not both at the same time.

The same goes for the integration of the second equation, you only get one integration constant that you can use to either fix $p_1(0)$ or $p_1(1)$, but it is improbable that you will get both to zero at the same time.

Usually a boundary value problem requires an order 2 ODE, then you would get 2 integration constants in every step that you can (try to) solve for the 2 boundary conditions.