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

% Solve driven pendulum equation phi''+alpha*phi'+sin(phi)=I order,
% constant step scheme), see Strogatz section 8.5, pages 265-273.

% define parameters:
% ------------------
alpha=1.2;
I=1.5;
delta_t=0.01;

%print (and save) every nprint steps:
nprint=1;
% integration time:
T=200.0;

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

% initial conditions:
x=zeros(2,NP);
x(:,1)=[0.0,0.1]';

% integrate using modified Euler scheme:
  for n=1:1:NP-1
    xp=x(:,n);
    for j=1:1:nprint
      x_tilde=xp+delta_t*rhs(xp,alpha,I);
      xp=xp+0.5*(rhs(xp,alpha,I)  ...
		 + rhs(x_tilde,alpha,I))*delta_t;
    end
    x(:,n+1)=xp;
  end


% start and end time for plotting:
b=1; %NP*3/4;
e=NP;

% plot simple time series:
time=delta_t*[1:1:NP];
plot(time(b:e),x(2,b:e))
xlabel('time')
ylabel('d phi/dt (rad.)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ff=rhs(x,alpha,I)
% return the rhs of the dynamical system:

ff=[x(2),I-sin(x(1))-alpha*x(2)]';
