連立常微分方程式をルンゲクッタ法で解く

以下のローレンツモデルを解く

clear; close all;

a=10; b=8/3; r=28;

f1 = @(t,x,y,z) -a*x+a*y;
f2 = @(t,x,y,z) -x*z+r*x-y;
f3 = @(t,x,y,z) x*y - b*z;

N = 2000;
t0 = 0;
tN = 40;
h = (tN-t0)/N;
t = linspace(t0,tN,N);

x = zeros(1,N);
y = zeros(1,N);
z = zeros(1,N);
x(1) = 1;
y(1) = 1;
z(1) = 1;

for k=2:N
  k1x = h*f1(t(k-1),x(k-1),y(k-1),z(k-1));
  k2x = h*f1(t(k-1)+h/2,x(k-1)+k1x/2,y(k-1)+k1x/2,z(k-1)+k1x/2);
  k3x = h*f1(t(k-1)+h/2,x(k-1)+k2x/2,y(k-1)+k2x/2,z(k-1)+k2x/2);
  k4x = h*f1(t(k-1)+h,x(k-1)+k3x,y(k-1)+k3x,z(k-1)+k3x);
  x(k) = x(k-1) + 1/6*(k1x+2*k2x+2*k3x+k4x);

  k1y = h*f2(t(k-1),x(k-1),y(k-1),z(k-1));
  k2y = h*f2(t(k-1)+h/2,x(k-1)+k1y/2,y(k-1)+k1y/2,z(k-1)+k1y/2);
  k3y = h*f2(t(k-1)+h/2,x(k-1)+k2y/2,y(k-1)+k2y/2,z(k-1)+k2y/2);
  k4y = h*f2(t(k-1)+h,x(k-1)+k3y,y(k-1)+k3y,z(k-1)+k3y);
  y(k) = y(k-1) + 1/6*(k1y+2*k2y+2*k3y+k4y);

  k1z = h*f3(t(k-1),x(k-1),y(k-1),z(k-1));
  k2z = h*f3(t(k-1)+h/2,x(k-1)+k1z/2,y(k-1)+k1z/2,z(k-1)+k1z/2);
  k3z = h*f3(t(k-1)+h/2,x(k-1)+k2z/2,y(k-1)+k2z/2,z(k-1)+k2z/2);
  k4z = h*f3(t(k-1)+h,x(k-1)+k3z,y(k-1)+k3z,z(k-1)+k3z);
  z(k) = z(k-1) + 1/6*(k1z+2*k2z+2*k3z+k4z);
end

figure();
plot(x,z);
xlabel('X');
ylabel('z');

figure();
plot(x,y);
xlabel('X');
ylabel('Y');


figure();
plot(t,x);
hold on;
plot(t,y);
plot(t,z);
legend('x','y','z');