%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% phase portrait for second order ODE
%
% mouse click in graphics window: start flow
%
% any button press at any position in window [ESC]: EXIT
%
% VHS Oct. 2003
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf
[y1,y2]=meshgrid(-4:0.5:8,-3:0.4:6);
v1=y2;
v2=-3*sin(y1)-0.2*y2;
quiver(y1,y2,v1,v2)
hold on
k = 0
while (k == 0)
    k = waitforbuttonpress;
    point = get(gca,'CurrentPoint');
    point = point(1,1:2);              % extract x and y
    [t,y]=ode23('damped',[0 10],point);
    plot(y(:,1),y(:,2),'r-');
end