clc;                                               % Clears the screen
clear all;

h=0.1;                                             % step size
x = 0:h:10;                                         % Calculates upto y(1)
y = zeros(1,length(x));
z = zeros(1,length(x));
y(1) = 1;                                          % initial condition
z(1) = 1;                                          % initial condition
% F_xy = @(t,r) 3.*exp(-t)-0.4*r;                  % change the function as you desire
F_xyz = @(x,y,z) z;                                  % change the function as you desire
G_xyz = @(x,y,z) -y ;

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

    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  % main equation
    z(i+1) = z(i) + (1/6)*(L_1+2*L_2+2*L_3+L_4)*h;  % main equation
    fprintf('%10d%+10.2f%+10.4f%+15.2fn', i, x(i+1), y(i+1), z(i+1));
end

plot(x, y, 'b', "linewidth", 2);

syms y(x)
Dy = diff(y);

ode = diff(y,x,2) ==  -y;
cond1 = y(1) == 1;
cond2 = Dy(1) == 1;
conds = [cond1 cond2];
ySol(x) = dsolve(ode,conds);
ySol = simplify(ySol);

hold on
fplot( ySol, [0, 10], "r", "linewidth", 2 )
hold off