%% Models the simple pendulum using the Improved Euler method. 
%% Usage: [t,th,w] = ImprovedEulerPendulum([ic1,ic2],tmin,tmax,N)
%% Input:  ic1 and ic2  - the intial values of theta (th) and w (which is th').
%%         tmin - initial time
%%         tmax - final time
%%         N - number of steps
%% Output: t - Nx1 vector of time values
%%         th - Nx1 vector of angle values
%%         w - Nx1 vector of angular velocity values
function [t,th,w] = ImprovedEulerPendulum(ic,tmin,tmax,N)
g = 9.81; L=1; % gravity and rod length
 
t = linspace(tmin,tmax,N)';
h = t(2) - t(1); % step size

th = zeros(size(t)); w = th;
th(1)=ic(1); w(1)=ic(2); % initial conditions

for j=1:N-1
   thdot1 = w(j); % slopes at previous points
   wdot1 = -g/L*sin(th(j));

   thdot2 = w(j) + h*wdot1; % slopes at Euler step
   wdot2 = -g/L*sin(th(j) + h*thdot1);
   
   th(j+1) = th(j) + h*(thdot1 + thdot2)/2;
   w(j+1) = w(j) + h*(wdot1 + wdot2)/2;
end

end