#{

% https://www.physicsforums.com/threads/problem-with-runge-kutta-script-file.1068182/

% all in one file:
% move Lornz function to the end  and add endfunction
% move rk4singlestep to the end   and add endfunction

#}

% the 3ed sript file :
 clear all ,close all , clc


% inputing the lorenz example parameters

sigma = 10;
beta = 8/3;
rho = 28;

% the initial condition of the lorenz function
y0 = [-8;8;27];

% computeing the trajectory
dt = 0.01;
tspan = 0:dt:4;
Y(:,1) = y0;
yk = y0;
for i=1:length(tspan)
time = i*dt;
ykplus1 = rk4singlestep(@(t,y)Lornz(t,y,sigma,beta,rho),dt,time,yk);
Y = [Y ykplus1];
yk = ykplus1;
end
plot3(Y(1,:),Y(2,:),Y(3,:),'b')


function dy = Lornz(t,y,sigma,beta,rho)
% y is the 3 dimentional state vector

% y(1) =x , y(2) = y , y(3) = z
dy = [
sigma*(y(2)-y(1));
y(1)*(rho-y(3))-y(2);
y(1)*y(2) - beta*y(3);
];
endfunction

% the second file rk4singlestep
% rang kuta function to solve the lorenz example or other function

function yout = rk4singlestep(fun,dt,tk,yk)
% fun is a function i should load into the memory from file or command
f1 = fun(tk,yk);
f2 = fun(tk+dt/2,yk+(dt/2)*f1);
f3 = fun(tk+dt/2, yk+(dt/2)*f2);
f4 = fun(tk+dt, yk + (dt)*f3);

yout = yk + (dt/6)*(f1+2*f2+2*f3+f4);

endfunction