function penHeun(N,t1)
%
% Use the Heun's method to integrate the
% equations of motion of a simple pendulum. 
% The values of the optional parameters N and t1
% determine the number of steps taken, and the 
% final time. The default values are 1000 and 30. 
%
% All rights reserved. Matt Kawski. April 1999.
% http://math.la.asu.edu/~kawski
%
if nargin < 2, t1=30; end
if nargin < 1, N=1000; end
%
% Define nominal values of the system parameters
% mass of the pendulum, length, grav accel, damping
m=1
L=9.81
g=9.81
beta=0
%
% Define the initial data
%
t0=0
theta0=(1/10)*pi
omega0=0
%
% Define the time interval on which to numerically 
% solve the DE. Specify the number of subintervals 
% to be used. 
% For efficiency reserve contiguous memory blocks 
% to store the values generated by the integration 
% algorithm. Define the first value in each array. 
%
dt=(t1-t0)/N
tt=[t0:dt:t1];
ttheta=zeros(size(tt));
oomega=zeros(size(tt));
ttheta(1)=theta0;
oomega(1)=omega0;
%
% Use Heun's method to solve the system of DEs.
%
for k=1:N
   a1=dt*(-beta/m/L*oomega(k)-g/L*sin(ttheta(k)));
   b1=dt*oomega(k);
   a2=dt*(-beta/m/L*(oomega(k)+a1)-g/L*(sin(ttheta(k)+b1)));
   b2=dt*(oomega(k)+a1);
   oomega(k+1)=oomega(k)+1/2*(a1+a2);
   ttheta(k+1)=ttheta(k)+1/2*(b1+b2);
end
%
% Plot the solutions.
%
plot(tt,oomega,tt,ttheta)
title('The simple pendulum via Heun''s method')
text(t1/10,1.2*max(ttheta),'Which curve is theta, which is omega?')
text(t1/10,1.1*max(ttheta),'Experiment with moderate N and large t1')