clear;  clc;

%   non linear pendulum, predictor-corrector

g=10;
L=0.5;
om=sqrt(g/L);
om2=om*om;

tstart=0;
tend=25;
dt=0.005;
y(1)=3.12;
t0=y(1);
y(2)=0;

k=0;
for time=tstart:dt:tend
    k=k+1;
    tt(k)=time;
    th(k)=y(1);
    true(k)=t0*cos(om*time);
    yprime=RHS(y,om2);
    ytilde = y + dt*yprime;
    yprimetilde=RHS(ytilde,om2);
    ytildetilde=y + dt*yprimetilde;
    y = 0.5*( ytilde + ytildetilde );
end

plot(tt,th,tt,true);

function dydt=RHS(y,om2)
    Theta=y(1);
    V=y(2);
    dydt(1)=V;
    dydt(2)=-om2*sin(Theta);
  % dydt(2)=-om2*Theta;
end