function ddp % Damped, driven pendulum
% Diff eq. is: d^2theta/dt^2 = -(g/l)*sin(theta)-b*dtheta/dt+FD*sin(omD*t)
%
% Rewrite as two, first order ODEs:
% dtheta/dt = y2
% dy1/dt = -(g/l)*sin(y1)-b*y2+FD*sin(omD*t)

set(0,'defaultaxesfontsize',20) % This sets the font size for the plots
close all
clc; clear

% Parameters are chosen to be in the "chaotic" regime
om2=1; % Natural frequency of the pendulum squared (g/l)
b = 0.5; % Friction strength
FD= 1.0; % Magnitude of the driving force
omD=2/3; % Driver frequency

y0 = [1.3, -1.6]; % Initial conditions

T = 2*pi/omD; % Period of the driver
N = 2405; % Number of periods
ts = 0:T:N*T; % Time span, stepping by T

% Make a theta vs. time plot
[t,y]=ode113(@f, [0 50], y0,[],om2,omD,FD,b);

plot(t,y(:,1),t,mod(y(:,1),2*pi),'.')
ylabel('\theta(t)')
xlabel('time')
title('Motion of a damped, driven pendulum')

figure % Open a second plot
[t,y]=ode113(@f, ts, y0,[],om2,omD,FD,b);

Nt = 5; % Number of "transient" steps to skip

% Make a poincare plot
plot(mod(y(Nt:end,1),2*pi),y(Nt:end,2),'.','markersize',1)
xlim([0,2*pi]);
ylim([-2, 1])
xlabel('\theta')
ylabel('\omega')
title('Poincare section for damped, driven pendulum')

function dydt=f(t,y,om2,omD,FD,b)
dydt = [y(2); -om2*sin(y(1))+FD*sin(omD*t)-b*y(2)]; % Compare to comments at top