%    dx/ds for particle paths
     function xdot=stagf(t,x)
     xdot=zeros(2,1);
     xdot(1)=x(1);
     xdot(2)=-x(2);
