Sorry, that should be Surfplot.  And the correct link (from the posting by jakubi in
<http://www.mapleprimes.com/questions/38398-surface-Plots>) is
<http://www.mapleprimes.com/view.aspx?sf=96307/280321/surfplot.zip>.  That is a zip archive containing files Surfplot.ind, Surfplot.lib and 8933_data2.mw.  Extract these to a directory, open the worksheet 8933_data2.mw, and follow the directions there.  Here's the plot I get:

Well, you won't get a symbolic sum of something like (-1)^n * n^(1/n), because there isn't a closed-form formula for such things (sorry, Maple doesn't know about the MRB constant).  In other cases where there is a formula, though, you might be able to get it using sum with the environment variable _EnvFormal set to true.  For example:

> _EnvFormal:= true:
   sum((-1)^n, n=1..infinity);

         -1/2

Essentially, Maple does use summability methods.  In this case, I think, Levin's U transform.

> evalf(Sum(n^(1/n)*(-1)^n, n=1..infinity));

                            -0.3121403575

> evalf(Sum((-1)^n,n=1..infinity));

                            -0.5000000000

Oops, that c[k] in A1 should have been c[j].  I just corrected it in my answer.

Oops, that c[k] in A1 should have been c[j].  I just corrected it in my answer.

@Kamel Algeria : "taken into account" maybe, but not by solve.  The error message is produced by assuming: Maple never gets around to calling solve.

@Kamel Algeria : "taken into account" maybe, but not by solve.  The error message is produced by assuming: Maple never gets around to calling solve.

@zhong chen : Preben plotted the function x(t), which is continuous.  It is the acceleration that is discontinuous, as I said.

@zhong chen : Preben plotted the function x(t), which is continuous.  It is the acceleration that is discontinuous, as I said.

Axel's code certainly should not give you a parse error or unterminated if.  What exactly did you try?

Axel's code certainly should not give you a parse error or unterminated if.  What exactly did you try?

It might help if you stated the whole question.  What are a, b, c, d and n?  What do you know about them?

It is possible to do this by finding a differential equation satisfied by the curves and plotting its solutions.
Note that F(x,y) = 0 satisfies the differential equation F_x + F_y dy/dx = 0 (where subscripts refer to partial derivatives).

First I get rid of the problem of complex values.

> F:=evalc(tan(x*sqrt(cp^2/3220^2-1))*(cp^2/3220^2-2)^2+4*tan(x*sqrt(cp^2/6450^2-1))
      *sqrt((cp^2/3220^2-1)*(cp^2/6450^2-1))) assuming cp > 3220, cp < 6450;

The result still contains signum and abs functions that should be simplified.


> Q:= op(op(indets(F,specfunc(anything,signum))));

Q := 1+1/431351361000000*cp^4-519709/4313513610000*cp^2

> signum(factor(Q)) assuming cp > 3220, cp < 6450;

                -1

 > F:= eval(F, {signum(Q)=-1, abs(Q) = -Q});

F := tan(1/3220*x*(cp^2-10368400)^(1/2))*(1/10368400*cp^2-2)^2-4*sinh(1/6450*x*(41602500-cp^2)^(1/2))*cosh(1/6450*x*(41602500-cp^2)^(1/2))/(1+sinh(1/6450*x*(41602500-cp^2)^(1/2))^2)*(-1-1/431351361000000*cp^4+519709/4313513610000*cp^2)^(1/2)

For our initial conditions, I'll find the x values at cp = 5000.

> ix:= select(type,[seq(fsolve(eval(F,cp=5000),x=i .. i+1),i=0..14)], float);

ix := [0., 1.250462379, 3.918863973, 6.564112774, 9.208734620, 11.85333446, 14.49793352]

Now for the DE

> de:= eval(diff(F,cp)+diff(F,x)*diff(X(cp),cp),x=X(cp));

I suspect the DE is stiff (because the curve comes so close to singularities), so I'll use a stiff solver.

> S1:= dsolve({de, X(5000)=0},numeric,stiff=true);

Now to solve with each initial condition.

> with(plots):
   for i from 1 to nops(ix) do 
      S1(initial=[5000, ix[i]]); PP[i]:= odeplot(S1,[X(cp),cp],cp=3220 .. 6450); 
   end do:

And to see the results:

> display([seq(PP[i],i=1..nops(ix))], view=[0..14, 3220 .. 6450]);

It is possible to do this by finding a differential equation satisfied by the curves and plotting its solutions.
Note that F(x,y) = 0 satisfies the differential equation F_x + F_y dy/dx = 0 (where subscripts refer to partial derivatives).

First I get rid of the problem of complex values.

> F:=evalc(tan(x*sqrt(cp^2/3220^2-1))*(cp^2/3220^2-2)^2+4*tan(x*sqrt(cp^2/6450^2-1))
      *sqrt((cp^2/3220^2-1)*(cp^2/6450^2-1))) assuming cp > 3220, cp < 6450;

The result still contains signum and abs functions that should be simplified.


> Q:= op(op(indets(F,specfunc(anything,signum))));

Q := 1+1/431351361000000*cp^4-519709/4313513610000*cp^2

> signum(factor(Q)) assuming cp > 3220, cp < 6450;

                -1

 > F:= eval(F, {signum(Q)=-1, abs(Q) = -Q});

F := tan(1/3220*x*(cp^2-10368400)^(1/2))*(1/10368400*cp^2-2)^2-4*sinh(1/6450*x*(41602500-cp^2)^(1/2))*cosh(1/6450*x*(41602500-cp^2)^(1/2))/(1+sinh(1/6450*x*(41602500-cp^2)^(1/2))^2)*(-1-1/431351361000000*cp^4+519709/4313513610000*cp^2)^(1/2)

For our initial conditions, I'll find the x values at cp = 5000.

> ix:= select(type,[seq(fsolve(eval(F,cp=5000),x=i .. i+1),i=0..14)], float);

ix := [0., 1.250462379, 3.918863973, 6.564112774, 9.208734620, 11.85333446, 14.49793352]

Now for the DE

> de:= eval(diff(F,cp)+diff(F,x)*diff(X(cp),cp),x=X(cp));

I suspect the DE is stiff (because the curve comes so close to singularities), so I'll use a stiff solver.

> S1:= dsolve({de, X(5000)=0},numeric,stiff=true);

Now to solve with each initial condition.

> with(plots):
   for i from 1 to nops(ix) do 
      S1(initial=[5000, ix[i]]); PP[i]:= odeplot(S1,[X(cp),cp],cp=3220 .. 6450); 
   end do:

And to see the results:

> display([seq(PP[i],i=1..nops(ix))], view=[0..14, 3220 .. 6450]);

@Christopher2222 : It is no safer to remove those 2944 draws than it is to remove any other set of 2944 draws.

Term to look up in your favourite probability text: "independent".