Troubles with ode45 solving two ODEs?

Question

I am solving two the first order ODEs:

$$p_0 p_0'=-\dfrac{32 \beta}{R^4}$$

$$(p_0p_1)'=-\dfrac{2-\sigma_v}{\sigma_v}\dfrac{8}{R}p_0'$$

I wrote following matlab code for solving it with ode45

function f= odefun_moja(Z,P)
P=zeros(2,1)
p0=P(1);
p1=P(2);
beta=1;
R=2;
sig=1;
%%%   EQUATIONS!!!
dp0dz=-32*beta/(p0*R^4);
dp1dz=(-((2-sig)/sig)*8*dp0dz/R-dp0dz*p1)/p0;
%%%   EQUATIONS!!!
f=[dpodz; dp1dz];

But after run with

function main_moja
z0=0;
zf=100;
zspan=[z0,zf];
y0=[1 0; 1 0];
options = odeset('RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));
[Z,P]=ode45(@odeset_moja,zspan,y0,options);

I got his message:

> Error using odearguments (line 92) ODEFUN_MOJA returns a vector of
> length 1, but the length of initial conditions vector is 4. The vector
> returned by ODEFUN_MOJA and the initial conditions vector must have
> the same number of elements.

When I try to solve first only the first equation, with code:

function f= odefun_moja(Z,P)
P=zeros(1,1)
p0=P(1);
beta=1;
R=2;
sig=1;
%%%   EQUATIONS!!!
dp0dz=-32*beta/(p0*R^4);
%%%   EQUATIONS!!!
f=[dp0dz];

Error message is this:

Error using feval
Error: File: odefun_moja.m Line: 15 Column: 1
Function definitions are not permitted in this context.

Error in odearguments (line 87)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

My questions:

1. Is that good approach for solving this system, start with the first equation, solve it for p0, and after that try to solve second equation? Or use both equations in code would be better?

2. I have next conditions $p_0|_{z=0}=p_{0i}$ and $p_0|_{z=1}=1$. According to literature it is necessary to find $p_0'|_{z=0}$ with shooting method, according to already mentioned $p_0|_{z=1}=1$. How to connect this two conditions and shoot $p_0'|_{z=0}$ for already known $p_0|_{z=1}=1$? Is that condition necessary, because this is first order equation, is there only one initial condition enough? (It will be the same process for second equation and $p_1$.)

Answer 1

You read the initialization thing wrong. Yes, you would have to initialize the return value

f=zeros(2,1)

if you then were to assign the components directly,

f(1) = -32*beta/(p0*R^4);                    %%% dp0dz
f(2) = (-((2-sig)/sig)*8*f(1)/R-f(1)*p1)/p0; %%% dp1dz

As you however directly construct the return array, no initialization of f is necessary. Note the typo in the first component.

---

You construct for y0 is strange, the solver expects a 2D vector. If you want to vary through different values of $p_1(0)$ you need to do that outside of the solver, with some kind loop over a list of those values.

---

It is unclear why you set zf = 100 if you want to integrate over the interval $[0,1]$.

---

In the given form with the given constants, your equations can be solved explicitly as $$ p_0(z)^2=p_0(0)^2-4z\\ p_0(z)p_1(z)=p_0(0)p_1(0)-4(p_0(z)-p_0(0)) $$ If the solution is to extend to $z=1$, the right side of the first equation needs to be positive there, thus $p_0(0)>2$. If you want $p_0(1)=1$ you need $p_0(0)=\sqrt5$. Now you can set either $p_1(0)$ or $p_1(1)$, the other of the two is then fixed by the chosen one.