Forced damped Pendulum

First, for a samplling of values of the parameter A, we plot simulations. We take samples of the solution once per period of the forcing using the events location feature of ode45.
b=1/2;
omega=2/3;
A_vector=[0.9 1.07 1.15 1.35 1.46 1.5];
behaviors=["Periodic orbit","Two-cycle","Chaos","Periodic Orbit","Four-cycle","Chaos"];
nA=length(A_vector);
y0=[0 1];
tSpan=[0 1000];
tTransient=200;
for k=1:nA
A=A_vector(k);
f=@(t,y)pendForceDamp(t,y,b,A,omega);
fEvent=@(t,y)pendulumEvent(t,y,omega);
options=odeset('Events',fEvent,'RelTol',1e-8,'AbsTol',1e-8);
[t,y,te,ye]=ode45(f,tSpan,y0,options);
% Throw away a transient
spot=find(t>tTransient,1,'first');
y=y(spot:end,:);
% Throw away a transient from the events time series
spot=find(te>tTransient,1,'first');
ye=ye(spot:end,:);
figure
plot(y(:,1),y(:,2),ye(:,1),ye(:,2),'o')
title(sprintf('$A=%0.2f$: %s',A),behaviors(k))
xlabel('$\theta$');ylabel('$\frac{d\theta}{dt}$')
figure
plot(mod(ye(:,1),2*pi),ye(:,2),'o')
axis([0 2*pi min(y(:,2)) max(y(:,2))])
title(sprintf('$A=%0.2f$: %s',A),behaviors(k))
xlabel('$\theta$');ylabel('$\frac{d\theta}{dt}$')
y0=y(end,:); % use the final solution as the new initial condtion
end