08-04-08.mw

 > 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)

 >