@blamm64 : Let me try to clarify. 

D can differentiate some procedures.  For example, here's a function of one variable which
implements 10 steps of the forward Euler method for solving the initial value problem
D(y)(t) = f(t, y(t)), y(0) = y0 (where f and y0 are global variables).

>  FE1:= proc(t)
       local h, i,  Y, T;
       Y:= y0; T:= 0;
       h:= t/10
       for i from 1 to 10 do
          Y:= Y + h*f(T,Y);
          T:= T + h;
       end do;
       Y
       end proc;

D can differentiate this, producing a new procedure that computes the derivative.  This is an exact symbolic calculation.

> D(FE1);

proc (t)
local T, Tt, Y, Yt, h, ht, i;
  Yt := 0;
  Y := y0;
  Tt := 0;
  T := 0;
  ht := 1/10;
  h := 1/10*t;
  for i to 10 do
    Yt := Yt+ht*f(T,Y)+h*(D[1](f)(T,Y)*Tt+D[2](f)(T,Y)*Yt);
    Y := Y+h*f(T,Y);
    Tt := Tt+ht;
    T := h+T
  end do;
  Yt
  end proc

If you change this to calculate the number of steps, D won't do it any more, but codegen[GRADIENT] will.

> FE2:= proc(t)
        local h, i, n, Y, T;
       Y:= y0; T:= 0;
       n:= ceil(100*t);
       h:= t/n;
       for i from 1 to n do
          Y:= Y + h*f(T,Y);
          T:= T + h;
       end do;
       Y
       end proc;

> D(FE2);

   D(FE2)

> codegen[GRADIENT](FE2);

  proc(t)
local T, Y, dT, dY, dh, dn, h, i, n;
    dY := 0;
    Y := y0;
    dT := 0;
    T := 0;
    dn := floor(1, 100*t);
    n := ceil(100*t);
    dh := 1/n - t*dn/n^2;
    h := t/n;
    for i to n do
        dY := h*D[1](f)(T, Y)*dT + (1 + h*D[2](f)(T, Y))*dY
             + f(T, Y)*dh;
        Y := Y + h*f(T, Y);
        dT := dT + dh;
        T := T + h
    end do;
    return dY
end proc

A separate issue is numeric differentiation, which evaluates a derivative numerically by evaluating the function at nearby points.  That's what evalf(D(f)(x)) will do if D can't differentiate f  symbolically.  This can be subject to severe roundoff problems.

As for documentation of eventstop, see the help page ?dsolve,numeric,Events. Yes, I think there is room for considerable improvement on this help page.

@blamm64 : Let me try to clarify. 

D can differentiate some procedures.  For example, here's a function of one variable which
implements 10 steps of the forward Euler method for solving the initial value problem
D(y)(t) = f(t, y(t)), y(0) = y0 (where f and y0 are global variables).

>  FE1:= proc(t)
       local h, i,  Y, T;
       Y:= y0; T:= 0;
       h:= t/10
       for i from 1 to 10 do
          Y:= Y + h*f(T,Y);
          T:= T + h;
       end do;
       Y
       end proc;

D can differentiate this, producing a new procedure that computes the derivative.  This is an exact symbolic calculation.

> D(FE1);

proc (t)
local T, Tt, Y, Yt, h, ht, i;
  Yt := 0;
  Y := y0;
  Tt := 0;
  T := 0;
  ht := 1/10;
  h := 1/10*t;
  for i to 10 do
    Yt := Yt+ht*f(T,Y)+h*(D[1](f)(T,Y)*Tt+D[2](f)(T,Y)*Yt);
    Y := Y+h*f(T,Y);
    Tt := Tt+ht;
    T := h+T
  end do;
  Yt
  end proc

If you change this to calculate the number of steps, D won't do it any more, but codegen[GRADIENT] will.

> FE2:= proc(t)
        local h, i, n, Y, T;
       Y:= y0; T:= 0;
       n:= ceil(100*t);
       h:= t/n;
       for i from 1 to n do
          Y:= Y + h*f(T,Y);
          T:= T + h;
       end do;
       Y
       end proc;

> D(FE2);

   D(FE2)

> codegen[GRADIENT](FE2);

  proc(t)
local T, Y, dT, dY, dh, dn, h, i, n;
    dY := 0;
    Y := y0;
    dT := 0;
    T := 0;
    dn := floor(1, 100*t);
    n := ceil(100*t);
    dh := 1/n - t*dn/n^2;
    h := t/n;
    for i to n do
        dY := h*D[1](f)(T, Y)*dT + (1 + h*D[2](f)(T, Y))*dY
             + f(T, Y)*dh;
        Y := Y + h*f(T, Y);
        dT := dT + dh;
        T := T + h
    end do;
    return dY
end proc

A separate issue is numeric differentiation, which evaluates a derivative numerically by evaluating the function at nearby points.  That's what evalf(D(f)(x)) will do if D can't differentiate f  symbolically.  This can be subject to severe roundoff problems.

As for documentation of eventstop, see the help page ?dsolve,numeric,Events. Yes, I think there is room for considerable improvement on this help page.

Another possibility is to enlarge the system you give to dsolve to include the quantity you're interested in: if you want to plot (or do something with) a certain quantity involving derivatives of your variables, add an equation

(that quantity) = q(t)

and dsolve will produce a solution that gives you the value of q(t).  Now if you do this for something involving a certain derivative of x, you'll need to provide initial values for lower derivatives of x.  If you find it difficult to calculate those: let dsolve provide them!  For example, consider the system

sys0:= {D(z)(t) = 3*y(t)-3, D(y)(t) = 2*x(t) - 2,
> D(x)(t) = 1/2*(z(t) - x(t)*y(t) + y(t) + x(t)),
> x(0)=1,y(0)=1,z(0)=1};

and say I want to plot q(t) = (D@@3)(x)(t)+(D@@2)(y)(t)+D(z)(t).  I'll need initial values D(x)(0), (D@@2)(x)(0), and D(y)(0).  To get the first derivatives:

> sys1:= sys0 union {D(x)(t) = xt(t), D(y)(t)=yt(t)};
   dsolve(sys1, numeric)(0);

To get the second derivative

> sys2:= sys1 union {(D@@2)(x)(t) = xtt(t), D(x)(0) = 1};
   dsolve(sys2, numeric)(0);

> sys3:= sys0 union {(D@@3)(x)(t)+(D@@2)(y)(t)+D(z)(t)=q(t), D(x)(0)=1, D(y)(0)=0,
(D@@2)(x)(0)=0};
    Sol:= dsolve(sys3, numeric);

> plots[odeplot](Sol, [t, q(t)], t = 0 .. 2);

Another possibility is to enlarge the system you give to dsolve to include the quantity you're interested in: if you want to plot (or do something with) a certain quantity involving derivatives of your variables, add an equation

(that quantity) = q(t)

and dsolve will produce a solution that gives you the value of q(t).  Now if you do this for something involving a certain derivative of x, you'll need to provide initial values for lower derivatives of x.  If you find it difficult to calculate those: let dsolve provide them!  For example, consider the system

sys0:= {D(z)(t) = 3*y(t)-3, D(y)(t) = 2*x(t) - 2,
> D(x)(t) = 1/2*(z(t) - x(t)*y(t) + y(t) + x(t)),
> x(0)=1,y(0)=1,z(0)=1};

and say I want to plot q(t) = (D@@3)(x)(t)+(D@@2)(y)(t)+D(z)(t).  I'll need initial values D(x)(0), (D@@2)(x)(0), and D(y)(0).  To get the first derivatives:

> sys1:= sys0 union {D(x)(t) = xt(t), D(y)(t)=yt(t)};
   dsolve(sys1, numeric)(0);

To get the second derivative

> sys2:= sys1 union {(D@@2)(x)(t) = xtt(t), D(x)(0) = 1};
   dsolve(sys2, numeric)(0);

> sys3:= sys0 union {(D@@3)(x)(t)+(D@@2)(y)(t)+D(z)(t)=q(t), D(x)(0)=1, D(y)(0)=0,
(D@@2)(x)(0)=0};
    Sol:= dsolve(sys3, numeric);

> plots[odeplot](Sol, [t, q(t)], t = 0 .. 2);

@blamm64 : D is the usual way to do symbolic differentiation of a procedure, but there's also codegen[GRADIENT].  Neither of these seem able to work on the procedures returned by dsolve.  For example, codegen[GRADIENT] returns errors such as

Error, (in intrep/statement) unable to translate Array(1..4, [...], datatype = anything)

I suspect there's way too much fancy stuff going on in these procedures to allow GRADIENT or D to work.  On the other hand, if you built a simple numerical procedure to implement one of the typical solution methods, D should work. 

@blamm64 : D is the usual way to do symbolic differentiation of a procedure, but there's also codegen[GRADIENT].  Neither of these seem able to work on the procedures returned by dsolve.  For example, codegen[GRADIENT] returns errors such as

Error, (in intrep/statement) unable to translate Array(1..4, [...], datatype = anything)

I suspect there's way too much fancy stuff going on in these procedures to allow GRADIENT or D to work.  On the other hand, if you built a simple numerical procedure to implement one of the typical solution methods, D should work. 

@hirnyk :

More generally, consider a function f(z) given by a power series

that has a positive radius of convergence, where m is a nonnegative integer.  Then

where

with the same radius of convergence.  Since g is continuous, we conclude that

as z -> 0.  Now apply this with z = 1/(2*n+1).

@hirnyk :

More generally, consider a function f(z) given by a power series

that has a positive radius of convergence, where m is a nonnegative integer.  Then

where

with the same radius of convergence.  Since g is continuous, we conclude that

as z -> 0.  Now apply this with z = 1/(2*n+1).

Actually it isn't the parser, rather it's a feature of Document Mode.  The document block containing the double for loop contains no output region.  When the worksheet is executed, e.g. with the !!! button, the code here is executed but its output is not shown.  If you put the cursor somewhere in the code and press Enter, an output region is produced, and you see the output (but even that won't happen if "Display output from Evaluate (Document Blocks)" is unchecked in Options,Display). 

I strongly recommend using Worksheet Mode rather than Document Mode, unless you really need to produce a Document.

Actually it isn't the parser, rather it's a feature of Document Mode.  The document block containing the double for loop contains no output region.  When the worksheet is executed, e.g. with the !!! button, the code here is executed but its output is not shown.  If you put the cursor somewhere in the code and press Enter, an output region is produced, and you see the output (but even that won't happen if "Display output from Evaluate (Document Blocks)" is unchecked in Options,Display). 

I strongly recommend using Worksheet Mode rather than Document Mode, unless you really need to produce a Document.

Your ball procedure needs its inputs to be numeric.  You're supplying it with 'f'(t) where t will be numeric, but because of the unevaluation quotes 'f'(t) will not.  It would work if you made sure that the arguments of ball were fully evaluated:

> ball := proc (x, y) plots[pointplot]([eval([x, y])], color = blue, symbol = solidcircle,
      symbolsize = 40) end proc;

Your ball procedure needs its inputs to be numeric.  You're supplying it with 'f'(t) where t will be numeric, but because of the unevaluation quotes 'f'(t) will not.  It would work if you made sure that the arguments of ball were fully evaluated:

> ball := proc (x, y) plots[pointplot]([eval([x, y])], color = blue, symbol = solidcircle,
      symbolsize = 40) end proc;

You might also look at the maptools package in the Application Center
http://www.maplesoft.com/applications/view.aspx?SID=3925

@blamm64 : yes, this is actually using numeric differentiation as D returns unevaluated on one of the solution procedures; accuracy could be a problem in some instances, but is probably OK here.

@blamm64 : yes, this is actually using numeric differentiation as D returns unevaluated on one of the solution procedures; accuracy could be a problem in some instances, but is probably OK here.