4th order runge-kutta method to solve two 1st order differential equations

%%4th order runge-kutta for Project 2
fy=@(t,y,z) z;
fz=@(t,y,z) (36000/(1500-40*t))-9.81-((0.1962*z^2)/(1500-40*t));

%boundary conditions
t(1) =0;
z(1) =0;
y(1)=0;
j=0;
h=1; %%step size
tfinal=10;
N=tfinal/h;
%Update loop

for j=1:N
    t(j+1)=t(j)+h;

    k1y=fy(t(j),y(j),z(j));
    k1z=fz(t(j),y(j),z(j));

    k2y=fy(t(j)+0.5*h , y(j)+0.5*h*k1y , z(j)+0.5*h*k1z);
    k2z=fz(t(j)+0.5*h,y(j)+0.5*h*k1y,z(j)+0.5*h*k1z);

    k3y=fy(t(j)+0.5*h,y(j)+0.5*h*k2y,z(j)+0.5*h*k2z);
    k3z=fz(t(j)+h/2,y(j)+0.5*h*k2y,z(j)+0.5*h*k2z);

    k4y=fy(t(j)+h,y(j)+h*k3y,z(j)+h*k3z);
    k4z=fz(t(j)+h,y(j)+h*k3y,z(j)+h*k3z);

    y(j+1)=y(j)+h/6*(k1y+2*k2y+2*k3y+k4y);

    z(j+1)=z(j)+h/6*(k1z+2*k2z+2*k3z+k4z);
end
fprintf('t = %6.4f, y = %6.4f\n', t, y);

so i need to display y values from t=0 to t=10 but the t output is not as expected. and y range shouldnt be a zero in it. please help

[SL changed the last line from a comment in the code to text in a text region.]

Answer 1

[Alan Stevens]

Like this?

%%4th order runge-kutta for Project 2
fy=@(t,y,z) z;
fz=@(t,y,z) (36000/(1500-40*t))-9.81-((0.1962*z^2)/(1500-40*t));

%boundary conditions

z(1) =0;
y(1)=0;

h=1; %%step size
tfinal=10;
N=tfinal/h;
%Update loop
t = 0:h:N*h;
for j=1:N
    

    k1y=fy(t(j),y(j),z(j));
    k1z=fz(t(j),y(j),z(j));

    k2y=fy(t(j)+0.5*h, y(j)+0.5*h*k1y, z(j)+0.5*h*k1z);
    k2z=fz(t(j)+0.5*h, y(j)+0.5*h*k1y, z(j)+0.5*h*k1z);

    k3y=fy(t(j)+0.5*h, y(j)+0.5*h*k2y, z(j)+0.5*h*k2z);
    k3z=fz(t(j)+0.5*h,  y(j)+0.5*h*k2y,  z(j)+0.5*h*k2z);

    k4y=fy(t(j)+h,  y(j)+h*k3y,  z(j)+h*k3z);
    k4z=fz(t(j)+h,  y(j)+h*k3y,  z(j)+h*k3z);

    y(j+1)=y(j)+h/6*(k1y+2*k2y+2*k3y+k4y);

    z(j+1)=z(j)+h/6*(k1z+2*k2z+2*k3z+k4z);
   
end
disp(['         t         y'])

         t         y

disp([t' y'])

         0         0
    1.0000    7.2008
    2.0000   29.2186
    3.0000   66.6541
    4.0000  120.0703
    5.0000  189.9819
    6.0000  276.8443
    7.0000  381.0419
    8.0000  502.8774
    9.0000  642.5612
   10.0000  800.2017

plot(t,y,'o-',t,z,'*-'),grid
xlabel('t'), ylabel('y & z')
legend('y','z')

[Jue]

yes thank you :)