I knew I had read about an event handler in Maple 13, but never had a chance to play with it until today. I don't claim to yet know how to make optimal use of this facility, but hope this post will encourage those with more knowledge to help me (and others) to improve our use of this intriguing new tool.

Here's the setup. I again use the two different versions of the last ODE in the system.

restart;
with( plots ):
G := 6.63*10^(-11)*10^(-33)*(1.49^3)(3600*(365*24))^12:
R := 2:
M := 3.6*10^37:
ic := r(0) = 7, v(0) = 1, omega(0) = 1, theta(0) = 1:
de1 := diff(r(t), t) = v(t):
de2 := diff(theta(t), t) = omega(t):
de3 := diff(omega(t), t) = -2*v(t)*omega(t)/r(t):
de4a := diff(v(t), t) = r(t)*omega(t)^2-G*M/r(t)^2; # if r(t)<10
de4b := diff(v(t), t) = r(t)*omega(t)^2-2*G*M*r(t)/R^3; # if r(t)>10

Note that I will create my own transition, at r(t)=10.

Since the IC starts with r(0)=7<10, the system initially follows the system with de4a:

s1 := dsolve([de1, de2, de3, de4a, ic], numeric, events=[[r(t)=10,halt]] ):
Warning, cannot evaluate the solution further right of .88446367, event #1 triggered a halt

Next, since the solution has r(t) increasing, the next stage uses de4b with the solution at this trigger point as the new initial conditions:

t1 := 0.88446367:
S1 := s1(t1):
ic1 := eval( S1[2..-1], S1[1] )[];
Warning, cannot evaluate the solution further right of .88446366, event #1 triggered a halt
omega(0.884463667843913392) = 0.490000011199886454,
r(0.884463667843913392) = 10.,
theta(0.884463667843913392) = 1.66681392288931306,
v(0.884463667843913392) = 5.06347526668841486
s2 := dsolve([de1, de2, de3, de4b, ic1], numeric, events=[[r(t)=10,halt]] ):
Warning, cannot evaluate the solution further right of 1.9665449, event #1 triggered a halt

The next stage in the solution uses de4a with new initial conditions:

t2 := 1.9665449:
S2 := s2(t2):
ic2 := eval( S2[2..-1], S2[1] )[];
Warning, cannot evaluate the solution further right of 1.9665448, event #1 triggered a halt
omega(1.96654486510690220) = 0.490000097951576841,
r(1.96654486510690220) = 10.0000000000000018,
theta(1.96654486510690220) = 2.11175299682023621,
v(1.96654486510690220) = -5.06347546127608882
s3 := dsolve([de1, de2, de3, de4a, ic2], numeric, events=[[r(t)=10,halt]] ):

I conclude by putting the three pieces together in a single plot:

P1 := odeplot( s1, [[t,r(t)],[t,omega(t)],[t,theta(t)],[t,v(t)]], t= 0..3, legend=[r,omega,theta,v] ):
P2 := odeplot( s2, [[t,r(t)],[t,omega(t)],[t,theta(t)],[t,v(t)]], t=t1..3, legend=[r,omega,theta,v] ):
P3 := odeplot( s3, [[t,r(t)],[t,omega(t)],[t,theta(t)],[t,v(t)]], t=t2..4, legend=[r,omega,theta,v] ):
display( [P1,P2,P3] );

The plot shows all four components of the solution, on each of the three time intervals.

I imagine there are better ways to implement this transition. Maybe even a way to do this automatically within the event handling system (see ?dsolve,numeric,Events ). I look forward to seeing how this could be done.

Doug

---------------------------------------------------------------------
Douglas B. Meade <><
Math, USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu
Phone: (803) 777-6183 URL: http://www.math.sc.edu