It's hard to comment without seeing the Maple program you're using.  For number-crunching, there are some very important ways to improve Maple's performance: in particular, use hardware floats, rtable-based data structures, and built-in or NAG routines as much as possible.  If you use DiscreteTransforms[FourierTransform] I think you'll get performance as good as, or maybe even better than, a Fortran or Pascal program you write yourself, but of course that's not really using "the Maple interpreter".

One problem is that the parameter c in animate is a float, but the index to plotter must be an integer from 1 to 15. 

There may also be a premature evaluation problem.

Try this:

> display([seq(pointplot3d(plotter[c],style=line,linestyle=solid,thickness=5),c=1..15)],insequence=true);

One problem is that the parameter c in animate is a float, but the index to plotter must be an integer from 1 to 15. 

There may also be a premature evaluation problem.

Try this:

> display([seq(pointplot3d(plotter[c],style=line,linestyle=solid,thickness=5),c=1..15)],insequence=true);

I don't see why you say the plot "should give the mean".  Your y simplifies to

ys := b/(1-x^2)^(1/2)/x/cosh(b*arccosh(1/x))^2*sinh(b*arccosh(1/x))

(note that it's arccosh, not acosh)

As x -> 1, the limit is b^2.  But I don't see why that should be in any sense a mean. Note in particular that the cosh(b*arccosh(1/x)) in the denominator will be 0 when  x = sec(Pi/(2*b)), which is in your x interval if b >= 12. Approaching this value of x, y will go to +infinity.

If your "simulated value" disagrees with this, perhaps there is something wrong with either the simulation or your equation.

The help pages ?plot,options and ?plot3d,options seem to indicate that only TIMES, COURIER, HELVETICA, and SYMBOL fonts are allowed.  But, at least in Maple 12 and 13  Standard Worksheet under Windows, arbitrary fonts actually do work.
 

> with(plots):
    textplot3d([2,3,4,"abcABC"],font=["Vivaldi',30],axes=box);

> textplot([seq([-i,1,cat("&#",1487+i,";")], i=1..27)], font=["Narkisim",30], axes=none),

The help pages ?plot,options and ?plot3d,options seem to indicate that only TIMES, COURIER, HELVETICA, and SYMBOL fonts are allowed.  But, at least in Maple 12 and 13  Standard Worksheet under Windows, arbitrary fonts actually do work.
 

> with(plots):
    textplot3d([2,3,4,"abcABC"],font=["Vivaldi',30],axes=box);

> textplot([seq([-i,1,cat("&#",1487+i,";")], i=1..27)], font=["Narkisim",30], axes=none),

A cute toy.  It rather reminds me of this:

www.theonion.com/content/video/apple_introduces_revolutionary

The use of dsolve(..., numeric) with parameters is described in the help page ?dsolve,numeric,interactive.  As Axel said, your solution will involve complex numbers, so I think you must use separate variables for the real and imaginary parts of x and z.

> eqs:= eval({eqn},{x(t)=xre(t)+I*xim(t),z(t)=zre(t)+I*zim(t)});
   eqs:=evalc({seq(seq(map(f,eqs[i]),i=1..2),f=[Re,Im])});
  inits:= {xre(0)=0,xim(0)=0,zre(0)=2.25,zim(0)=0,D(xre)(0)=0,D(xim)(0)=0,
      D(zre)(0)=0.5,D(zim)(0)=0};
  S:= dsolve(eqs union inits,numeric,parameters=[g], maxfun=0);

> S(parameters=[1]);

> S(1);

[t = 1., xim(t) = -.987164142564391224, diff(xim(t),t) = -1.99881117720197166, xre(t) = -4.27344649929832254, diff(xre(t),t) = -8.65780250921212336, zim(t) = 4.27587766787401336, diff(zim(t),t) = 8.65780239210555180, zre(t) = 1.26033041347707498, diff(zre(t),t) = -1.99893730354143706]

> plots[odeplot](S,[[t,xre(t)],[t,xim(t)],[t,zre(t)],[t,zim(t)]],t=0..1,colour=[red,green,blue,black]);

The use of dsolve(..., numeric) with parameters is described in the help page ?dsolve,numeric,interactive.  As Axel said, your solution will involve complex numbers, so I think you must use separate variables for the real and imaginary parts of x and z.

> eqs:= eval({eqn},{x(t)=xre(t)+I*xim(t),z(t)=zre(t)+I*zim(t)});
   eqs:=evalc({seq(seq(map(f,eqs[i]),i=1..2),f=[Re,Im])});
  inits:= {xre(0)=0,xim(0)=0,zre(0)=2.25,zim(0)=0,D(xre)(0)=0,D(xim)(0)=0,
      D(zre)(0)=0.5,D(zim)(0)=0};
  S:= dsolve(eqs union inits,numeric,parameters=[g], maxfun=0);

> S(parameters=[1]);

> S(1);

[t = 1., xim(t) = -.987164142564391224, diff(xim(t),t) = -1.99881117720197166, xre(t) = -4.27344649929832254, diff(xre(t),t) = -8.65780250921212336, zim(t) = 4.27587766787401336, diff(zim(t),t) = 8.65780239210555180, zre(t) = 1.26033041347707498, diff(zre(t),t) = -1.99893730354143706]

> plots[odeplot](S,[[t,xre(t)],[t,xim(t)],[t,zre(t)],[t,zim(t)]],t=0..1,colour=[red,green,blue,black]);

The second derivative of f at x is expressed either as diff(f(x),x$2) or as (D@@2)(f)(x).

Do you mean the limit of y as x -> 1?

> limit(b*sech(b*arccosh(1/x))*tanh(b*arccosh(1/x))/(sqrt((1/x)-1)*sqrt(1+(1/x))*(x^2)),x=1);

b^2

Or for more detail,

> map(simplify, series(b*sech(b*arccosh(1/x))*tanh(b*arccosh(1/x))/(sqrt((1/x)-1)*sqrt(1+(1/x))*(x^2)),
     x=1));

b^2+5*b^2*(-1+b^2)*(x-1)*(1/3)+(1/30)*b^2*(64-125*b^2+61*b^4)*(x-1)^2+(1/126)*b^2*(-312+889*b^2-854*b^4+277*b^6)*(x-1)^3+(1/22680)*b^2*(-229000*b^2+324093*b^4+50521*b^8+62136-207750*b^6)*(x-1)^4+O((x-1)^5)

Actually I'd call that two periods.

Since _EnvAllSolutions:= true causes solve(cos(x)=1) to return all the solutions (2*Pi*_Z2), why does it not do so for the inequality?

Actually I'd call that two periods.

Since _EnvAllSolutions:= true causes solve(cos(x)=1) to return all the solutions (2*Pi*_Z2), why does it not do so for the inequality?

It seems to me that odeplot requires a single result of dsolve(...,numeric), so that would mean only one IC (unless you make several copies of the system with different IC's into a single larger system).  However, you could use display to  combine several odeplot results.

As for critical points: the points I had there were the ones with x' = 0, but you could modify it:

> ts:= [seq(fsolve(subs(sols[i],x(t)*(1-0.5*x(t)-0.5*y(t))),0..5), i=1..2)];
  pts:= [seq(subs(sols[i],[x(t),y(t)])(ts[i]),i=1..2)];
  P2:= pointplot(pts,symbol=solidcircle,symbolsize=16,colour=black): 
  P3:= arrow([pts[1]+[1.1,1.1],pts[2]+[1.1,1.1]],[<-1,-1>,<-1,-1>]):
  display([P1,P2,P3],scaling=constrained,view=[0..4,-1..4]);

It seems to me that odeplot requires a single result of dsolve(...,numeric), so that would mean only one IC (unless you make several copies of the system with different IC's into a single larger system).  However, you could use display to  combine several odeplot results.

As for critical points: the points I had there were the ones with x' = 0, but you could modify it:

> ts:= [seq(fsolve(subs(sols[i],x(t)*(1-0.5*x(t)-0.5*y(t))),0..5), i=1..2)];
  pts:= [seq(subs(sols[i],[x(t),y(t)])(ts[i]),i=1..2)];
  P2:= pointplot(pts,symbol=solidcircle,symbolsize=16,colour=black): 
  P3:= arrow([pts[1]+[1.1,1.1],pts[2]+[1.1,1.1]],[<-1,-1>,<-1,-1>]):
  display([P1,P2,P3],scaling=constrained,view=[0..4,-1..4]);