h = 0.03125;      % Stepsize
tspan = [0 20];   % Time interval
y0 = [5 5 5];     % Initial conditions
a = 10;           
b = 8/3;     
r = 28;
[tl, yl] = rk4sys(@lorenz,tspan,y0,h,a,b,r); % RK4 solution
% Time-series plot using RK4 method
figure(1)
plot(tl,yl(1:end,1))
legend('x = y = z = 5')
title('RK4 time plot')
xlabel("Time")
ylabel("Quantity")

function [tp,yp] = rk4sys(dydt,tspan,y0,h,varargin)
% Reference: Chapra, S. C. (2018). Applied Numerical Methods
% with MATLAB® for Engineers and Scientists Fourth Edition, pp.602-603
% rk4sys: fourth-order Runge-Kutta for a system of ODEs
% [t,y] = rk4sys(dydt,tspan,y0,h,p1,p2,...): integrates
% a system of ODEs with fourth-order RK method
% input:
% dydt = name of the M-file that evaluates the ODEs
% tspan = [ti, tf]; initial and final times with output
% generated at interval of h, or
% = [t0 t1 ... tf]; specific times where solution output
% y0 = initial values of dependent variables
% h = step size
% p1,p2,... = additional parameters used by dydt
% output:
% tp = vector of independent variable
% yp = vector of solution for dependent variables
if nargin < 4,error('at least 4 input arguments required'), end
if any(diff(tspan)<= 0),error('tspan not ascending order'), end
n = length(tspan);
ti = tspan(1);tf = tspan(n);
if n == 2
t = (ti:h:tf)'; n = length(t);
if t(n)<tf
t(n + 1) = tf;
n = n + 1;
end
else
t = tspan;
end
tt = ti; y(1,:) = y0;
np = 1; tp(np) = tt; yp(np,:) = y(1,:);
i = 1;
while(1)
tend = t(np + 1);
hh = t(np + 1) - t(np);
if hh > h,hh = h;end
while(1)
if tt+hh > tend,hh = tend-tt;end
k1 = dydt(tt,y(i,:),varargin{:})';
ymid = y(i,:) + k1*hh/2;
k2 = dydt(tt + hh/2,ymid,varargin{:})';
ymid = y(i,:) + k2*hh/2;
k3 = dydt(tt + hh/2,ymid,varargin{:})';
yend = y(i,:) + k3*hh;
k4 = dydt(tt + hh,yend,varargin{:})';
phi = (k1 + 2*(k2 + k3) + k4)/6;
y(i + 1,:) = y(i,:) + phi*hh;
tt = tt + hh;
i = i + 1;
if tt >= tend,break,end
end
np = np + 1; tp(np) = tt; yp(np,:) = y(i,:);
if tt >= tf,break,end
end

function yl = lorenz(t,y,a,b,r)
yl = [-a*y(1)-a*y(2);
    r*y(1)-y(2)-y(1)*y(3);
    -b*y(3)+y(1)*y(2)];
end

function yl = lorenz(t,y,a,b,r)
yl = [-a*y(1)+a*y(2);
    r*y(1)-y(2)-y(1)*y(3);
    -b*y(3)+y(1)*y(2)];
end