MATLAB: Differential eguations: modify the program

differential eguations

Differential eguations
dy1/dt=-a*y1+y2*y3
dy2/dt=-b(y2-y3)
dy3/dt=-a*y2+ y3-3*y2*y3
a=4 b=8 c=30
y1(0)=0
y2(0)=0
y3(0)=10-10 =0.0000000001
MATLAB ODE45_main:
tspan=[0,100];
y0=[0; 0;10^(-10)];
[t,y]=ode45('ODE45_fun',tspan,y0);
data=[t,y];
save ODE45_data.txt data –ascii
plot3(y(:,1), y(:,2), y(:,3))
grid on
MATLAB ODE45_fun:
function dy=ODE45_fun(t,y)
a=4; b=18; c=30;
dy(1)=-a*y(1)+y(2)*y(3);
dy(2)=-b*(y(2)-2*y(3));
dy(3)=c*y(2)-y(3)-3*y(2)*y(1);
dy=[dy(1);dy(2);dy(3)];
The program can work without any error. If the formula for y4:
y4=-3*y2*y3-4*y1*y3+8*y1*y3
How can I modify the program? And the curve of y4 ?

Best Answer

tspan=[0,100];
y0=[0; 0;10^(-10);0];
M=[1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 0];
options=odeset('Mass',M);
[t,y]=ode45(@ODE45_fun,tspan,y0,options);
function dy=ODE45_fun(t,y)
a=4; b=18; c=30;
dy(1)=-a*y(1)+y(2)*y(3);
dy(2)=-b*(y(2)-2*y(3));
dy(3)=c*y(2)-y(3)-3*y(2)*y(1);
dy(4)=y(4)-(-3*y(2)*y(3)-4*y(1)*y(3)+8*y(1)*y(3));
dy=[dy(1);dy(2);dy(3);dy(4)];
Best wishes
Torsten.