% FP Forced pendulum.

% Francisco J. Beron-Vera, 2006/15/05

clear, clf

options=odeset('RelTol',1e-6,'AbsTol',1e-6);

% Parameters.
a=1; c=0; f=0; omega=0;
%a=1; c=0; f=.5; omega=1; T=2*pi/omega;
%a=1; c=.05; f=5; omega=1; T=2*pi/omega;

% Initial conditions.
y0=[-5:-3 .5:.5:2 3:5]';  
x0=repmat(2*pi,size(y0));
z0=[x0';y0']; 
z0=z0(:);

% Integration.
if omega==0
   [t z]=ode45(@fp2eq,[0 200*2*pi],z0(:),...
      options,...
      a,c,f,omega);
   x=z(:,1:2:end-1);
   y=z(:,2:2:end);
else
   [t z]=ode45(@fp2eq,[0 50*T],z0(:),...
      options,...
      a,c,f,omega);
   x=z(:,1:2:end-1);
   y=z(:,2:2:end);
   [ti z]=ode45(@fp2eq,0:T:500*T,z0(:),...
      options,...
      a,c,f,omega);
   xi=z(:,1:2:end-1);
   yi=z(:,2:2:end);
end

% Ploting...
subplot(221)
plot(t,y)
if omega==0
   axis([0 10*2*pi -5 5])
else
  axis([0 10*T -5 5]) 
end
xlabel('$t$',...
   'Interpreter','latex','FontSize',14)
ylabel('$y$',...
   'Interpreter','latex','FontSize',14)
subplot(222)
plot(mod(x,4*pi),y,'.','MarkerS',1)
axis([pi 3*pi -5 5])
set(gca,'XTick',[pi 2*pi 3*pi],'XTickL',[-3.14 0 3.14])
if omega==0
   xlabel('$x$ mod $2\pi$',...
      'Interpreter','latex','FontSize',14)
end
title('Phase Portrait',...
   'Interpreter','latex','FontSize',14)
if omega~=0
   subplot(224)
   plot(mod(xi,4*pi),yi,'.','MarkerS',1)
   axis([pi 3*pi -5 5])
   set(gca,...
      'XTick',[pi 2*pi 3*pi],...
      'XTickL',[0 3.14 6.28])
   xlabel('$x$ mod $2\pi$',...
      'Interpreter','latex','FontSize',14)
   title('Poincar\''e Section',...
      'Interpreter','latex','FontSize',14)
end