Runge Kutta 5th order help, error message "Index exceeds matrix dimensions."

Hello I have written the following code for a second order differential equation, split into two first order differential equations using the 5th order Runge Kutta method but keep getting the error message "Index exceeds matrix dimensions.". Could someone please explain how to resolve this error? My original equations are:

dt/dx=z

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

Start of code:

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)
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]

You don't have any x(i+1)=etc line. Since you don't set the next value of x, when you get to the second loop iteration the x(2) reference fails.

[Jeviah Forrest]

thank you that was an oversight on my part, I've been reviewing the results and am unsure as to what else is wrong the the code. All values of z and t are 0 after the first step which should not be the case, any idea as to why that is?

[James Tursa]

You don't show your current code, but I am assuming you have added this line?

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

For the t and z variables, I would have expected the updates to be added to the current value, not multipled by the current value. Also, it appears you have a confusion with the variable h. At the top of your code and in your derivative equations it is used as a parameter in the calculations, but in your update code it appears you are using it as a stepsize (I assume the h here is because that is the standard variable name used for stepsize in textbooks etc). 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;

Note: It is a lot easier to read code when you put in spacing.

[Jeviah Forrest]

Thank you very much the above steps solved my problem, much appreciated and I will take note of spacing my code.