function doall=pendulum
clear all; close all;
% Eli Tziperman: example for course in nonlinear dynamics and
% chaos.

% Solve simple pendulum problem numerically using the "modified Euler scheme"
% (a second order, constant step scheme), see Strogatz page 33.

% the equation:
% d^2\theta/dt^2=(g/L)*\sin(\theta) - gamma*(d\theta/dt) + f*cos(omega_f*t)

%  First, translate to dynamical system of first order ODEs:
%  dx/dt = y
%  dy/dt = -(g/L)*sin(x) - gamma*y + f*cos(omega_f*z)
%  dz/dt = 1

% define parameters:
% ------------------
% pendulum parameters: gravity and length:
g=9.8;
L=9.8;
% Linear frequency:
omega_L=sqrt(g/L);
% Strength of friction:
gamma=0.1;
% f is forcing amplitude:
% Seems that increasing f to 0.8 or such will get chaotic solution;
% 0.5 gives a periodic solution;  0.0 shows decaying oscillations at ...
%     natural frequency which looks nice in phase space:
f=0.0;
% Forcing frequency:
omega_f=0.5345678*sqrt(g/L);
% time step: needs to be small enough for a given problem:
delta_t=0.01*2*pi/(max(omega_f,omega_L));
% Print (and save) every nprint steps:
nprint=1;
% Integration time:
T=1500.0;

% Number of steps:
N=floor(T/delta_t);
NP=floor(N/nprint);

% Initial conditions:
x=zeros(3,NP);
n=1;
x(:,n)=[pi/8.0,0,0]';

% Integrate:
for n=1:1:NP
    xp=x(:,n);
    for j=1:1:nprint
         x_tilde=xp+delta_t*rhs(xp,g,L,gamma,omega_f,f);
         xp=xp+0.5*(rhs(xp,g,L,gamma,omega_f,f)  ...
      +                rhs(x_tilde,g,L,gamma,omega_f,f))*delta_t;
    end
    x(:,n+1)=xp;
end

% Plot simple time series:
plot(x(3,:),x(1,:))
xlabel('time')
ylabel('theta (rad.)')

delay=floor(1.0/(delta_t*sqrt(g/L)));
b=floor(NP/2);
e=NP-3*delay;

% Plot phase space:
%  figure
%  plot3(x(1,b:e),x(2,b:e),x(3,b:e))
%  %plot(x(1,b:e),x(2,b:e))
%  xlabel('theta')
%  ylabel('d(theta)/dt')

% Plot reconstructed phase space using delay coordinates:
figure
plot3(x(2,b:e),x(2,b+delay:e+delay),x(2,b+2*delay:e+2*delay))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ff=rhs(x,g,L,gamma,omega_f,f)
% return the rhs of the dynamical system:

ff=[x(2),-(g/L)*sin(x(1)) - gamma*x(2) + f*cos(omega_f*x(3)),1.0]';