# Numerical solution of ODEs system using ODE45

2 views (last 30 days)
Maria Cuellar on 16 Feb 2020
Edited: James Tursa on 17 Feb 2020 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.
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);

James Tursa on 16 Feb 2020
Edited: James Tursa on 17 Feb 2020
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');