% Surface of Section Plot of damped, forced Duffing equation
%
% d^2 y/dt^2 = y-y^3-delta*dy/dt + A*cos(omega t)

% Use inline function
clear

delta=0.25;
gamma=0.30;
om=1.0;
fduff = @(t,y) [y(2); y(1)-y(1)^3-delta*y(2)+gamma*cos(om*t)];

[t,x]=ode45(fduff,[0:2*pi/om:5000/om],[1,0]);
plot(x(:,1),x(:,2),'.','MarkerSize',1)
axis([-2 2 -2 2])
