An attempt at question 2 of the xkcd Velociraptor Math problem (mentioned on this blog post). The parameters and events facilities of dsolve,numeric are handy here.

P[x],P[y]:=6*t*cos(phi),6*t*sin(phi):

ic1:=R1[x](0)=0,R1[y](0)=20*sqrt(3)*(1/3):
ic2:=R2[x](0)=10,R2[y](0)=-10*sqrt(3)*(1/3):

S1:=piecewise(t<10*(1/4),4*t,10):
de1x:=(D(R1[x]))(t)
      = S1*(P[x]-R1[x](t))/sqrt((P[x]-R1[x](t))^2+(P[y]-R1[y](t))^2):
de1y:=(D(R1[y]))(t)
      = S1*(P[y]-R1[y](t))/sqrt((P[x]-R1[x](t))^2+(P[y]-R1[y](t))^2):

S2:=4*t:
de2x:=(D(R2[x]))(t)
      = S2*(P[x]-R2[x](t))/sqrt((P[x]-R2[x](t))^2+(P[y]-R2[y](t))^2):
de2y:=(D(R2[y]))(t)
      = S2*(P[y]-R2[y](t))/sqrt((P[x]-R2[x](t))^2+(P[y]-R2[y](t))^2):

evt1:=(P[x]-R1[x](t))^2+(P[y]-R1[y](t))^2-0.1e-9:
evt2:=(P[x]-R2[x](t))^2+(P[y]-R2[y](t))^2-0.1e-9:

numsolver:=dsolve({ic1,ic2,de1x,de1y,de2x,de2y},'numeric',
            'parameters'=[phi],'events'=[[evt1,'halt'],[evt2,'halt']]):

interface('warnlevel'=0):

objective:=proc(phi_start)
global numsolver;
  numsolver('parameters'=[phi=phi_start]);
  numsolver(1000);
  op([1,2],numsolver('last'));
end proc:

sol:=Optimization:-Maximize(objective,-Pi/2..Pi/2):

tm,best:=sol[1],sol[2][1]; # best time, and angle in radians

convert(best,'units','radians','degrees'); # best angle in degrees

with(plots): with(plottools):
p1,p2:=map(rhs,[ic1]),map(rhs,[ic2]):
display(odeplot(numsolver,eval([
                          [R1[x](t),R1[y](t),'color'='red'],
                          [R2[x](t),R2[y](t),'color'='blue'],
                          [P[x],P[y],'color'='green','linestyle'='dash']
                          ],phi=best),
                0..tm,'legend'=["Raptor 1","Raptor 2","person"]),
        line(p1,p2),line(p2,[0,p2[2]]),'scaling'='constrained');

The results (in Maple 13.02) of 3.1 seconds and 32.8 degrees seem to agree roughly with a few others (here, and in the discussion thread pages).

Suggested improvements are welcome.