Numerical solution of ODEs system using ODE45

First I did the variable sustition as follow:

x=x(1)

x'=x(2)

x''=-x(1)-x(4)-3*x(2)^2+x(5)^3+6*x(6)^2+2*t

y=x(4)

y'=x(5)

y''=x(6)

y'''=-x(6)-x(2)-exp(-x(1))-t

The code is:

ecu1=@(t,x)[x(2);-x(1)-x(4)-3*x(2)^2+x(5)^3+6*x(6)^2+2*t;x(4);x(5);x(6);-x(6)-x(2)-exp(-x(1))-t]
rvt=0:0.001:5;
ci=[2 -4 -2 7 6];
[t,x]=ode45(ecu1,rvt,ci)
plot(t,x)

Index exceeds array bounds.

Error in Estudio1>@(t,x)[x(2);-x(1)-x(4)-3*x(2)^2+x(5)^3+6*x(6)^2+2*t;x(4);x(5);x(6);-x(6)-x(2)-exp(-x(1))-t]

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

Thanks for your help

Answer 1

[James Tursa]

You've got a 2nd order DE and a 3rd order DE, so the order of your system is 2+3=5, not 6. So your state vector should have 5 elements, not 6 elements. Your state definitions should be:

x=x(1)

x'=x(2)

y=x(3)

y'=x(4)

y''=x(5)

Which means your derivatives should be

dx(1) = x' = x(2)

dx(2) = x'' = -x(1) - x(3) - (3*x(2))^2 + x(4)^3 + 6*x(5) + 2*t

dx(3) = y' = x(4)

dx(4) = y'' = x(5)

dx(5) = y''' = -x(5) - x(2) - exp(-x(1)) - t

So your derivative function should be this instead of what you have

ecu1=@(t,x)[x(2);-x(1)-x(3)-(3*x(2))^2+x(4)^3+6*x(5)+2*t;x(4);x(5);-x(5)-x(2)-exp(-x(1))-t]

Also, the initial conditions seem to start at t=1, so that would imply that the first time of your rvt variable should be 1, not 0.

Finally, you probably don't want to plot the entire x variable since that contains all of the x', y', and y'' stuff as well. Maybe something like this instead to get x and y plotted:

plot(t,x(:,1),t,x(:,3));
xlabel('t');
legend('x','y');