Here's an "exact" (and much faster) way to do the cone projection.

ConeProj2:= proc(r,h,A,B,C)
 uses LinearAlgebra;
 local V,P,M,d;
 V:= evalf(GramSchmidt([<A,B,C>,<1,0,0>,<0,1,0>,<0,0,1>], conjugate=false,normalized=true));
 P:= Matrix(<V[2]^%T,V[3]^%T>);
 M:= r * P[1..2,1..2];
 d:= Norm(M^(-1).P.<0,0,h>,2);
 if d <= 1 then evalf(Pi*abs(Determinant(M)))
 else evalf((sqrt(d^2-1) + Pi - arcsin(sqrt(d^2-1)/d))*abs(Determinant(M)));
 end if
end proc;

Here's an "exact" (and much faster) way to do the cone projection.

ConeProj2:= proc(r,h,A,B,C)
 uses LinearAlgebra;
 local V,P,M,d;
 V:= evalf(GramSchmidt([<A,B,C>,<1,0,0>,<0,1,0>,<0,0,1>], conjugate=false,normalized=true));
 P:= Matrix(<V[2]^%T,V[3]^%T>);
 M:= r * P[1..2,1..2];
 d:= Norm(M^(-1).P.<0,0,h>,2);
 if d <= 1 then evalf(Pi*abs(Determinant(M)))
 else evalf((sqrt(d^2-1) + Pi - arcsin(sqrt(d^2-1)/d))*abs(Determinant(M)));
 end if
end proc;

@hirnyk : Ah, it seems to need another evalf.

ConeProj:= proc(r,h,A,B,C)
    uses LinearAlgebra;
    local V,verts;
    V:= evalf(GramSchmidt([<A,B,C>,<1,0,0>,<0,1,0>,<0,0,1>], conjugate=false,normalized=true));
    verts:= map(v -> [DotProduct(v,V[2],conjugate=false),
                               DotProduct(v,V[3],conjugate=false)],
       evalf([<0,0,h>,seq(<r*cos(2*Pi*t),r*sin(2*Pi*t),0>,t=0 .. 0.999,.001)]));
    simplex[convexhull](verts, output=area);
  end proc;

@hirnyk : Ah, it seems to need another evalf.

ConeProj:= proc(r,h,A,B,C)
    uses LinearAlgebra;
    local V,verts;
    V:= evalf(GramSchmidt([<A,B,C>,<1,0,0>,<0,1,0>,<0,0,1>], conjugate=false,normalized=true));
    verts:= map(v -> [DotProduct(v,V[2],conjugate=false),
                               DotProduct(v,V[3],conjugate=false)],
       evalf([<0,0,h>,seq(<r*cos(2*Pi*t),r*sin(2*Pi*t),0>,t=0 .. 0.999,.001)]));
    simplex[convexhull](verts, output=area);
  end proc;

It's really hard to guess what's going wrong without seeing your equations and the command you're using.  Please post them (or upload a worksheet).

@tatan2501 :
<http://books.google.ca/books?id=BsSwAAAAIAAJ>

@tatan2501 :
<http://books.google.ca/books?id=BsSwAAAAIAAJ>

@Christopher2222 : unfortunately I didn't keep the code for many of those animations; those I do have are not pretty, as they were just intended as "one-off" jobs for my own use only.  Here's a worksheet for one of the "predator-prey" animations.

ppreyn.mws

Oh, and how could I leave out the all-time classic? W. Feller, "An Introduction to Probability Theory and Its Applications", volumes I and II.

Oh, and how could I leave out the all-time classic? W. Feller, "An Introduction to Probability Theory and Its Applications", volumes I and II.

@tatan2501 : Yes, they were made with Maple (in Maple versions going back as far as Maple V Release 3).   In some cases I think I may have made individual .gif files with Maple and put them together into an animation using ImageMagick.

You really need some context to make sense of typical Maple mathematical animations.  FWIW, here's an old web page of mine on the Double Pendulum, including two animations made with Maple.

<http://www.math.ubc.ca/~israel/m215/doublpen/doublpen.html>

And, from the same course, some more pages with Maple animations,

The Snowplow Problem: <http://www.math.ubc.ca/~israel/m215/plows/plows.html>

Forced Vibrations: <http://www.math.ubc.ca/~israel/m215/forced/forced.html>

Phase Portraits of Linear Systems: <http://www.math.ubc.ca/~israel/m215/linphase/linphase.html>

Predator and Prey: <http://www.math.ubc.ca/~israel/m215/predprey/predprey.html>

See the help page ?Optimization,General,Methods.

If the method option is not provided, the method is chosen by NLPSolve according to the following rules.  If the optimization problem is univariate and unconstrained except for finite bounds, quadratic interpolation is used.  If the problem is unconstrained and the gradient of the objective function is available, the PCG method is used.  Otherwise, the SQP method is used.

See the help page ?Optimization,General,Methods.

If the method option is not provided, the method is chosen by NLPSolve according to the following rules.  If the optimization problem is univariate and unconstrained except for finite bounds, quadratic interpolation is used.  If the problem is unconstrained and the gradient of the objective function is available, the PCG method is used.  Otherwise, the SQP method is used.

Applying Preben's corrections, and adding an additional boundary condition f(150) = 1, another error occurs:

> ics := f(0) = 0, D(f)(150) = 0, D(f)(0) = 0, f(150) = 1;
   res := dsolve({ics, ode1}, numeric);

Error, (in dsolve/numeric/bvp) system is singular at left endpoint, use midpoint method instead

Maple is right: the system is indeed singular at the left endpoint because of the boundary condition f(0) = 0 and the fact that the highest-order derivative is multiplied by f(eta) in the differential equation.  But the error message is slightly strange, because there is no "midpoint" method to use here.  What should be used, I think, are method=bvp[midrich] or method = bvp[middefer].  However, both of those also produce an error:

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

Perhaps a suitable approximate solution could be provided with the approxsoln option, or a continuation method could be used.  Or perhaps there simply isn't a solution to this boundary value problem.