>	xphug:= [ diff(theta(t),t) = ( v(t)^2 - cos(theta(t))) / v(t),         diff(v(t),t)     = -sin(theta(t)) - R*v(t)^2 ,         diff(x(t),t)     = v(t)*cos(theta(t)),         diff(y(t),t)     = v(t)*sin(theta(t))]; with(DETools):with(plots):

 

(1)

 

>

bisector := proc( f, low, high, epsilon)   local mid, lo, hi;   lo := evalf(low);   hi := evalf(high);   if (f(lo)>0) then error("f not neg at lo=",lo, "f=",f(lo)); fi;   if (f(hi)<0) then error("f not pos at hi=",hi, "f=",f(hi)); fi;   mid := evalf((lo+hi)/2);   while ( abs(f(mid))>epsilon) do      if ( f(mid) > 0) then         hi:= mid;      else         lo:= mid;      fi;      mid := (lo+hi)/2; print( [lo, hi], mid, "f(mid)=", f(mid));   od;   return(mid); end:

 

>

R:=0.3;

 

(2)

 

>	sol:=dsolve({op(xphug), theta(0)=0, v(0)=2, x(0)=0, y(0)=1},                numeric,  range=0..10,output=procedurelist);

 

(3)

 

>

sol(4);

 

(4)

 

>

sol(5);

 

(5)

 

>	saul:=dsolve({op(xphug), theta(0)=0, v(0)=2, x(0)=0, y(0)=1},                numeric,  range=0..10, output=listprocedure);

 

(6)

 

>	yt:=saul[5];

 

(7)

 

>

yt(5);

 

(8)

 

>

sol(4);

 

(9)

 

>

saul(4);

 

(10)

 

>	sol[5](4);

 

(11)

 

>	saul[5](4);

 

(12)

 

>	DEplot(xphug, [ theta(t), v(t), x(t), y(t) ], t=0..8,

 

>	[seq([theta(0)=(Pi*(i/10)), v(0)=2, x(0)=0, y(0)=1 ], i=-3..2)],

 

>	theta=-Pi/2..3*Pi, v=0..3, x=-1..6, y=-1..3,         scene=[x,y], title="path of glider",

 

>	linecolor=[cyan,blue,violet,red,yellow,green]);

 

>

 

>	moses:=dsolve({op(xphug), theta(0)=-Pi/5, v(0)=2, x(0)=0, y(0)=1},                numeric,  range=0..10, output=listprocedure);

 

(13)

 

>	moses[5](8);

 

(14)

 

>	moses[5](6);

 

(15)

 

>	rhs(roses=red);

 

(16)

 

>	rhs(moses[5](5));

 

(17)

 

>	plot(rhs(moses[5](t)),t=4..6);

 

>	rhs(moses[5](5.22));

 

(18)

 

>	rhs(moses[5](5.21));

 

(19)

 

>	bisector(t->rhs(moses[5](t)), 5.22, 5.21, .000001);

 

 

 

 

 

 

 

 

 

 

	(20)

 

>	moses[4](5.215869142);

 

(21)

 

>	xphug:= [ diff(theta(t),t) = piecewise(y(t)>0,( v(t)^2 - cos(theta(t))) / v(t),0),         diff(v(t),t)     = piecewise(y(t)>0,-sin(theta(t)) - R*v(t)^2 ,0),         diff(x(t),t)     = piecewise(y(t)>0, v(t)*cos(theta(t)), 0),         diff(y(t),t)     = piecewise(y(t)>0, v(t)*sin(theta(t)), 0)];

 

(22)

 

>	crash:=dsolve({op(xphug), theta(0)=-Pi/5, v(0)=2, x(0)=0, y(0)=1},                numeric,  range=0..10, output=listprocedure);

 

(23)

 

>	crash[5](8);

 

(24)

 

>

crash(8);

 

(25)

 

>	splat:=Theta->dsolve({op(xphug), theta(0)=Theta, v(0)=2, x(0)=0, y(0)=1},                numeric,  range=0..10, output=listprocedure)(10);

 

(26)

 

>	splat(-Pi/5);

 

(27)

 

>	splat(-Pi/5+.01);

 

(28)

 

>