Incorrect results using Runge Kutta 5th order method

Hello, I have written the following code for a second order differential equation, split into two first order differential equations, the code itself seems to work however after the first step all results for z and t are 0 which is incorrect, is there any errors in my code/working?

Original equation: d^2t/dx^2 = (h*p/k*ac)*(t-ta)

Split equations: dt/dx = z

dz/dx = (h*p/k*ac)*(t-ta)

5th order Runge Kutta:

ac=7.854*10^-7 %cross sectional area

k=110; % Conductive heat transfer coefficient

ta=25; % Room temperature

h=20; %Convective heat transfer coefficient

p=0.04063; %Perimeter of cylinder

s=0.001 %step size

xfinal=0.02; % final distance

N=xfinal/s;

x(1)=0; %initial values

t(1)=50;

z(1)=0;

dtdx=@(x,t,z) z;

dzdx=@(x,t,z) (h*p/k*ac)*(t-ta);

for i=1:N

k1= dtdx(x(i),t(i),z(i));

l1=dzdx(x(i),t(i),z(i));

k2=dtdx(x(i)+0.25*s,t(i)+0.25*k1*s,z(i)+0.25*l1*s)

l2=dzdx(x(i)+0.25*s,t(i)+0.25*k1*s,z(i)+0.25*l1*s)

k3=dtdx(x(i)+0.25*s,t(i)+0.125*k1*s+0.125*k2*s,z(i)+0.125*l1*s+0.125*l2*s)

l3=dzdx(x(i)+0.25*s,t(i)+0.125*k1*s+0.125*k2*s,z(i)+0.125*l1*s+0.125*l2*s)

k4=dtdx(x(i)+0.5*s,t(i)-0.5*k2*s+k3*s,z(i)-0.5*l2*s+l3*s)

l4=dzdx(x(i)+0.5*s,t(i)-0.5*k2*s+k3*s,z(i)-0.5*l2*s+l3*s)

k5=dtdx(x(i)+0.75*s,t(i)+0.1875*k1*s+0.5625*k4*s,z(i)+0.1875*l1*s+0.5625*l4*s)

l5=dzdx(x(i)+0.75*s,t(i)+0.1875*k1*s+0.5625*k4*s,z(i)+0.1875*l1*s+0.5625*l4*s)

k6=dtdx(x(i)+s,t(i)-0.4286*k1*s+0.2857*k2*s+1.7143*k3*s-1.7143*k4*s+1.1429*k5*s,z(i)-0.4286*l1*s+0.2857*l2*s+1.7143*l3*s-1.7143*l4*s+1.1429*l5*s)

l6=dzdx(x(i)+s,t(i)-0.4286*k1*s+0.2857*k2*s+1.7143*k3*s-1.7143*k4*s+1.1429*k5*s,z(i)-0.4286*l1*s+0.2857*l2*s+1.7143*l3*s-1.7143*l4*s+1.1429*l5*s)

x(i+1)=x(i)+s;

t(i+1)=t(i)*0.0111*(7*k1+32*k3+12*k4+32*k5+7*k6)*h

z(i+1)=z(i)*0.0111*(7*l1+32*l3+12*l4+32*l5+7*l6)*h

end

Answer 1

[James Tursa]

I have already answered this on your other thread, but will re-post here. The 'h' variable in your t and z updates is supposed to be stepsize, but you are confusing it with your heat parameter that is used in your derivative equations. So you probably need to use 's' here. Also, I would expect your updates to be added to the current values, not multiplied by your current values. So maybe something like this is what you need:

t(i+1) = t(i) + 0.0111*( 7*k1 + 32*k3 + 12*k4 + 32*k5 + 7*k6 )*s;
z(i+1) = z(i) + 0.0111*( 7*l1 + 32*l3 + 12*l4 + 32*l5 + 7*l6 )*s;