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

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

Answer 1

[James Tursa]

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)

Answer 2

[Manoj Kumar Ghosh]

clc

clear all

n=10000;

x=linspace(1,3,n+1);

y0=[0;2/3];

h=(x(end)-x(1))/n;

f=@(x,y,z) ([z;((5*z)/(3*x)-(5*y)/(3*x.^2))]) %took y'=z

for i=1:n

k1=h*f(x(i),y0(1),y0(2));

k2=h*f(x(i)+h/2,(y0(1)+k1(1)/2),(y0(2)+k1(2)/2));

k3=h*f(x(i)+h/2,(y0(1)+k2(1)/2),(y0(2)+k2(2)/2));

k4=h*f(x(i)+h,(y0(1)+k3(1)),(y0(2)+k3(2)));

y0=y0+(1/6)*(k1+2*k2+2*k3+k4);

end

disp(y0(1))