function penRK4(N,t1)
%
% Use a 4th order Runge-Kutta 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 a 4th order Runge Kutta 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/2)-g/L*(sin(ttheta(k)+b1/2)));
   b2=dt*(oomega(k)+a1/2);
   a3=dt*(-beta/m/L*(oomega(k)+a2/2)-g/L*(sin(ttheta(k)+b2/2)));
   b3=dt*(oomega(k)+a2/2);
   a4=dt*(-beta/m/L*(oomega(k)+a3/2)-g/L*(sin(ttheta(k)+b3/2)));
   b4=dt*(oomega(k)+a3);
   oomega(k+1)=oomega(k)+1/6*(a1+2*a2+2*a3+a4);
   ttheta(k+1)=ttheta(k)+1/6*(b1+2*b2+2*b3+b4);
end
%
% Plot the solutions.
%
plot(tt,oomega,tt,ttheta)
title('The simple pendulum via 4th order Runge Kutta')
text(t1/10,1.1*max(ttheta),'Note, even with RK4 the amplitudes still grow, albeit very slowly')