Runge-Kutta method for ODE system of 3 equations

Question

I am trying to solve the Thomas system using the 4th order Runge-Kutta method, my code so far looks like this:

clc;                                             
clear all;
b=0.18;
h=0.001;                                             % step size
k = 0:h:1;  
x = zeros(1,length(k));                       
y = zeros(1,length(k)); 
z = zeros(1,length(k)); 
x(1) = 1;
y(1) = 0;                                          % initial condition
z(1) = 1;                                          % initial condition
F_xyz = @(k,x,y,z) sin(y)-b*x;                                  
G_xyz = @(k,x,y,z) sin(z)-b*y;
H_xyz = @(k,x,y,z) sin(x)-b*z;

for i=1:(length(x)-1)                              % calculation loop
k_1 = F_xyz(k(i),x(i),y(i),z(i));
L_1 = G_xyz(k(i),x(i),y(i),z(i));
Q_1 = H_xyz(k(i),x(i),y(i),z(i));
k_2 = F_xyz(k(i)+0.5*h,x(i)+0.5*h*k_1,y(i)+0.5*h*L_1,z(i)+0.5*h*Q_1);
L_2 = G_xyz(k(i)+0.5*h,x(i)+0.5*h*k_1,y(i)+0.5*h*L_1,z(i)+0.5*h*Q_1);
Q_2 = H_xyz(k(i)+0.5*h,x(i)+0.5*h*k_1,y(i)+0.5*h*L_1,z(i)+0.5*h*Q_1);
k_3 = F_xyz((k(i)+0.5*h),(x(i)+0.5*h*k_2),(y(i)+0.5*h*L_2),  (z(i)+0.5*h*Q_2));
L_3 = G_xyz((k(i)+0.5*h),(x(i)+0.5*h*k_2),(y(i)+0.5*h*L_2),(z(i)+0.5*h*Q_2));
Q_3 = H_xyz((k(i)+0.5*h),(x(i)+0.5*h*k_2),(y(i)+0.5*h*L_2),(z(i)+0.5*h*Q_2));
k_4 = F_xyz((k(i)+h),(x(i)+k_3*h),(y(i)+L_3*h),(z(i)+Q_3*h));      
L_4 = G_xyz((k(i)+h),(x(i)+k_3*h),(y(i)+L_3*h),(z(i)+Q_3*h));
Q_4 = H_xyz((k(i)+h),(x(i)+k_3*h),(y(i)+L_3*h),(z(i)+Q_3*h));

z(i+1) = z(i) + (1/6)*(Q_1+2*Q_2+2*Q_3+Q_4)*h;
x(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  % main equation
y(i+1) = z(i) + (1/6)*(L_1+2*L_2+2*L_3+L_4)*h;  % main equation

end

But this code does not provide the correct plot (I am using plot3 to plot x,y,z) as it should show the Thomas attractor for b=0.18. Any advice?

Answer 1

In the last three lines the updating of $x$ and $y$ is not correct (you update $x$ using $y$ and $y$ using $z$).