function poincare(T,a,b)
%This function M-file computes the Poincare map for the equation x'=f(t,x) on the
%interval [a,b], where f(t,x) is the function defined in the function M-file field.m, assumed
%T-periodic in t. Executing poincare produces a graph of the Poincare map on [a,b],
%and of the identity map.
x0=linspace(a,b);
for k=1:100
     p(k)=poinc(x0(k),T);
end
%This for loop is needed since the input x0 in ode45 cannot be a vector.
figure
plot(x0,[x0;p])
xlabel('x(0)')
ylabel('x(T)')
title('Poincare map')
grid on
shg
%The function y=poinc(x0,T) solves the IVP for x'=f(t,x) on [0,T],
%with IC x(0)=x0. y is the scalar x(T).
function y=poinc(x0,T)
sol=ode45(@field,[0,T],x0);
y=deval(sol,T);
%written and tested by Alex Freire, 9/4/05. Should work with Matlab 6.0 or later versions.