Index exeeded error for ode45 function

hi all,

I have generated a ODE matrix (M)as shown in the code.

ri=0.53;

syms w(r)
Nrad(r)=-(( (r^-2) - 1)./((ri^-2)+1));
Ntheta(r)=(( (r^-2) + 1)./((ri^-2)-1));

[V] = odeToVectorField(diff(w, 4) + (2/r)*diff(w,3) - ((1/r^2) + 6*Nrad(r) )*diff(w,2) + ((((1/r^2)-6*Ntheta(r) ))/r)*diff(w,1)  == 0 )
M = matlabFunction(V,'vars', {'r','Y'})
sol = ode45(M,[0.53 0.9],[0 0.53]);
fplot(@(x)deval(sol,x,1), [0.53, 0.9])

The ode45 is generating a error as follows:

what is causing this?

Index exceeds the number of array elements (2).
Error in
symengine>@(r,Y)[Y(2);Y(3);Y(4);-(1.0./r.^2.*3.157935826372082e-1-1.315793582637208).*Y(3)-(Y(4).*2.0)./r+((1.0./r.^2.*1.343763037129746+2.343763037129746).*Y(2))./r]

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

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

  Error in ODE_AnnularDisc_Nagarjuna (line 17)
sol = ode45(M,[0.53 0.9],[0 0.53]);

Answer 1

[Stephan]

This is an ode of order 4 - so odeToVectorField will reformulate it to a system of 4 ode with order 1:

>> pretty(V)
/                           Y                            \
|                            2                           |
|                                                        |
|                           Y                            |
|                            3                           |
|                                                        |
|                           Y                            |
|                            4                           |
|                                                        |
| /   3221    5618 \                                     |
| | ------- + ---- | Y                                   |
| |       2   2397 |  2   2 Y                            |
| \ 2397 r         /         4   /   4045     16854 \    |
| --------------------- - ---- - | -------- - ----- | Y  |
|           r               r    |        2   12809 |  3 |
\                                \ 12809 r          /    /

Thus you have to tell ode45 a number of 4 initial points - one for every equation. Using some fantasy values would give:

ri=0.53;

syms w(r)
Nrad(r)=-(( (r^-2) - 1)./((ri^-2)+1));
Ntheta(r)=(( (r^-2) + 1)./((ri^-2)-1));

[V] = odeToVectorField(diff(w, 4) + (2/r)*diff(w,3) - ((1/r^2) + 6*Nrad(r) )*diff(w,2) + ((((1/r^2)-6*Ntheta(r) ))/r)*diff(w,1)  == 0 )
M = matlabFunction(V,'vars', {'r','Y'})
sol = ode45(M,[1 10],[0 0 0.53 0.9]);
fplot(@(x)deval(sol,x,1), [1, 10])