MATLAB: Runge-Kutta 4th Order on a 2nd Order ODE

calculus

I am dealing with the following 2nd order ODE:
With initial conditions of: , , .
and should be 3.24
I have found the system as:
For the life of me I cannot see to get the below code to produce the correct output… is it the code or my calculus?
clear all;
fy=@(x,y,z) 6*x*z-5*z
fz=@(x,y,z) z;
x(1)=0;
z(1)=2/3;
y(1)=0;
h=0.5;
xfinal=3;
N=ceil((xfinal-x(1))/h);
for j=1:N
    x(j+1)=x(j)+h;
    k1y=fy(x(j),y(j),z(j));
    k1z=fz(x(j),y(j),z(j));
    k2y=fy(x(j)+h/2,y(j)+h/2*k1y,z(j)+h/2*k1z);
    k2z=fz(x(j)+h/2,y(j)+h/2*k1y,z(j)+h/2*k1z);
    k3y=fy(x(j)+h/2,y(j)+h/2*k2y,z(j)+h/2*k2z);
    k3z=fz(x(j)+h/2,y(j)+h/2*k2y,z(j)+h/2*k2z);
    k4y=fy(x(j)+h,y(j)+h*k3y,z(j)+h*k3z);
    k4z=fz(x(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(y(N));
figure;
plot(x,y,'LineWidth',2);
xlabel('X');
ylabel('Y');

Best Answer

All of your derivatives should be with respect to the independent variable x. In the context of your problem, dy/dz does not make any sense. So you have two variables y and z, and two derivatives dy/dx and dz/dx. So your two 1st order equations are:
dy/dx = y' = z
dz/dx = d(y')/dx = y'' = (5*x*y' - 5*y) / (3*x^2) = (5*x*z - 5*y) / (3*x^2)
The dz/dx expression is obtained by setting your differential equation to 0 (you didn't specify this but I assumed it to be the case) and then solving for y'' and substituting in z for y'. So your function handles should be:
fy=@(x,y,z) z
fz=@(x,y,z) (5*x*z - 5*y) / (3*x^2)
Related Question