% Solve nonlinear pendulum motion using midpoint method

g = 9.8; % Gravity
L = 1; % Length

theta0 = 45; % Initial angle in degrees
omega0 = 0;  % Initial angular rotation rate (degrees/sec)

Tmax = 10; % Final time
dt = 0.15;  % Time step

% Calculate number of steps
N = ceil(Tmax/dt) + 1;
t = 0:dt:((N-1)*dt);
v = zeros(2,N);
f = zeros(2,1);

% Set initial condition
v(:,1) = pi/180*[theta0; omega0];

% Use forward Euler for first step
f(1) = v(2,1); 
f(2) = -g/L*sin(v(1,1));
  
% Use forward Euler
v(:,2) = v(:,1) + dt*f;

% Perform calculation
for n = 2:N-1,
  
  f(1) = v(2,n); 
  f(2) = -g/L*sin(v(1,n));
  
  % Use midpoint method
  v(:,n+1) = v(:,n-1) + 2*dt*f;
  
end

plot(t,v(1,:)*180/pi);