% Solutions of pendulum

theta0=10*pi/180/2;
k=sin(theta0/2);
m=k^2;
K=ellipke(m);

g=9.8;
L=1.0;
omega=sqrt(g/L);

t=0:0.01:10;

figure(1)
plot(t,2*asin(k*ellipj(omega*t+K,k)),'k',t,theta0*sin(omega*t+pi/2),'k--')
legend('Nonlinear','Linear','location','southwest')
axis([0,10,-1.5*theta0,1.5*theta0])

figure(2)
[SN,CN,DN] = ellipj(omega*t+K,m);
theta=2*asin(k*SN);
plot(t,theta);
hold

%Numerical Solution
per=2*pi*sqrt(L/g);
a=0;
%b=5*per;
b=10;
N=1000;
h=(b-a)/N;

% Initial Condition
t2(1)=0;
theta2(1)=theta0;
omega2(1)=0;

% Euler-Cromer Method
for i=2:N+1
    omega2(i)=omega2(i-1)-g/L*h*sin(theta2(i-1));
    theta2(i)=theta2(i-1)+h*omega2(i);
    t2(i)=t2(i-1)+h;
end

plot(t2+h,theta2,'r');
hold off