MATLAB: Having issues with RK4 on small angle pendulum

Question

I am trying to implemnt and compaire some numerical solvers, but having issues when it comes to the RK4 as its not behaving as i would expect.

I would expect it to be stable and follow the analytic soultion, but it becomes unstable. And it cant figure out whats wrong in my rk4 code.

equations to be solved:
omega_dot = -(g/l)*theta
theta_dot=omega;
clf
clear all
%Declearing stuff
l = 1;
m = 1;
g = 9.81;
theta(1) = 0.2;
omega(1) = 0;
delta_t = 0.01;
time = 10; 
steps = time/delta_t;
%Analytic answer
t = [0:delta_t:time];
theta_a = theta(1)*cos(sqrt(g/l)*t);
%Euler solver
for i = 1:steps
   E_E(i)= (1/2)*m*l^2*omega(i)^2+m*g*l*(1-cos(theta(i))); %#ok<*SAGROW>

   omega(i+1) = omega(i)-(g/l)*theta(i)*delta_t; 
   theta(i+1) = theta(i)+omega(i)*delta_t;
   
end
%Last step of total Energy

E_E (i+1)= (1/2)*m*l^2*omega(i+1)^2+m*g*l*(1-cos(theta(i+1)));
figure(1)
hold on;
% plot(t,theta)

plot(t,theta_a,'.')
figure(2)
hold on;
plot(t,E_E)
%Euler-Cromer
for i = 1:steps
   E_EC (i)= (1/2)*m*l^2*omega(i)^2+m*g*l*(1-cos(theta(i))); %#ok<*SAGROW>
   omega(i+1) = omega(i)-(g/l)*theta(i)*delta_t; 
   theta(i+1) = theta(i)+omega(i+1)*delta_t;
   
end
%Last step of total Energy
E_EC (i+1)= (1/2)*m*l^2*omega(i+1)^2+m*g*l*(1-cos(theta(i+1)));
% 
% figure(1)
% hold on;
% plot(t,theta)
% figure(2)
% plot(t,E_EC)
for i = 1:steps
    
    k1 = -(g/l)*theta(i);
    k2 = -(g/l)*(theta(i)+delta_t*0.5*k1);
    k3 = -(g/l)*(theta(i)+delta_t*0.5*k2);
    k4 = -(g/l)*(theta(i)+delta_t*k3);
    omega(i+1) = omega(i)+((k1+2*k2+2*k3+k4)/6)*delta_t;
    
    k1 = omega(i);
    k2 = omega(i)+delta_t*0.5*k1;
    k3 = omega(i)+delta_t*0.5*k2;
    k4 = omega(i)+delta_t*k3;
    theta(i+1) = theta(i)+((k1+2*k2+2*k3+k4)/6)*delta_t; 
      
    
end
figure(1)
hold on;
plot(t,theta,'x')

Answer 1

Please use button fro code inserting