function peneuler(N,t1)
%
% Use the plain Euler'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 a plain Euler's method to solve the system of DEs.
%
for k=1:N
   oomega(k+1)=oomega(k)+dt*...
                  (-beta/m/L*oomega(k)-g/L*sin(ttheta(k)));
   ttheta(k+1)=ttheta(k)+dt*oomega(k);
end
%
% Plot the solutions.
%
plot(tt,oomega,tt,ttheta)
title('The simple pendulum via Euler''s method')
text(t1/10,6/5*max(ttheta),'Which curve is theta, which is omega?')
text(t1/10,max(ttheta),'Are the increasing amplitudes reasonable?')
text(t1/10,4/5*max(ttheta),'Where do they come from? Try increasing N')