Having issues with RK4 on small angle pendulum

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

[darova]

Please use button fro code inserting

[Truls Kjøsnes Olsen]

Yeah sorry about the code!

And thanks alot :)

[darova]

I think each step of RK4 should be connected

for i = 1:steps-1
    k1o = -(g/l)*theta(i);
    k1t = omega(i);
    k2o = -(g/l)*(theta(i)+delta_t*0.5*k1t);
    k2t = omega(i) + delta_t*0.5*k2o;
    k3o = -(g/l)*(theta(i)+delta_t*0.5*k2t);
    k3t = omega(i) + delta_t*0.5*k3o;
    
    k4o = -(g/l)*(theta(i)+delta_t*k3t);
    k4t = omega(i) + delta_t*k4o;
    
    omega(i+1) = omega(i)+((k1o+2*k2o+2*k3o+k4o)/6)*delta_t;
    theta(i+1) = theta(i)+((k1t+2*k2t+2*k3t+k4t)/6)*delta_t;  
end

Dimensions must agree