%    dx/ds for streamlines
     function xdot=stream(s,x)
     global tstream
     xdot=zeros(2,1);
     xdot(1)=x(1)/(1+tstream);
     xdot(2)=x(2);
