>

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),         diff(blackbox(t),t) = piecewise(y(t)>0, 1, 0)]: with(DETools):with(plots):

 

>

R:=0.3;

 

(1)

 

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

 

>	crash(10);

 

(2)

 

>	crash(5.21589191691218712);

 

(3)

 

>	crash(5.2);

 

(4)

 

>	f:= proc(Theta)  local  crash;  crash:=dsolve({op(xphug), theta(0)=Theta, v(0)=2, x(0)=0, y(0)=1, blackbox(0)=0},                numeric,  range=0..10, output=listprocedure):  return(rhs(crash(10)[5])); end:

 

>

 

>

f(-Pi/5);

 

(5)

 

>	plot(f,-Pi/4..Pi/4);

 

>

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:

 

>

f(-.61);

 

(6)

 

>

f(-.65);

 

(7)

 

>

f(.50);

 

(8)

 

>

f(.48);

 

(9)

 

>	bisector(t->f(t)-5, -.61, -.65, .00001);

 

Error, (in bisector) f not neg at lo=, -.61, f=, 0.47661354e-1

 

>	bisector(t->f(t)-5, -.65, -.61, .00001);

 

 

 

 

 

 

 

 

 

 

 

 

 

	(10)

 

>	bisector(t->f(t)-5, .50, .48, .00001);

 

 

 

 

 

 

 

 

 

 

 

	(11)

 

>