The display command is in the plots package.  So you could use plots[display] instead of just display.

The option insequence=true will give you an animation.  Leave it out if you want the different solutions all on one graph.

The display command is in the plots package.  So you could use plots[display] instead of just display.

The option insequence=true will give you an animation.  Leave it out if you want the different solutions all on one graph.

1) Your X and G are tables.  ScatterPlot needs Vectors or lists.  I'm surprised you don't get an "invalid input" error message.
You might try

> ScatterPlot(convert(X,list), convert(G,list));

2) To restrict a plot P to a certain window, x from a to b and y from c to d (where a,b,c,d are real constants):

> plots[display](P, view=[a..b, c..d]); 

I think Dirk was talking about parentheses around the C:

which results when

plot(x^2,x=0..5,labels=[x,theta [°C]]);

is copied into 2D input.  This seems to interpret the second label as theta*[`°`(C).

What I would do is type it as theta*[`°`*C] in Maple notation input.

I think Dirk was talking about parentheses around the C:

which results when

plot(x^2,x=0..5,labels=[x,theta [°C]]);

is copied into 2D input.  This seems to interpret the second label as theta*[`°`(C).

What I would do is type it as theta*[`°`*C] in Maple notation input.

I'm guessing that .9999999996 is really supposed to be 1,  1.380952380 is 29/21, 1.166666667*10^(-14) is 7/(10*b).

Then your expression L[1] is

> f := cos(7*Pi*x/(10*b))*cos(3*Pi*x/(10*b)) - 29/21 * sin(7*Pi*x/(10*b))*sin(3*Pi*x/(10*b));

> combine(f);

-4*cos(2*Pi*x/(5*b))*(1/21)+25*cos(Pi*x/b)*(1/21)

Now if w = exp(I*Pi*x/(5*b)), this can be written as

> g := -2/21*w^2 - 2/21*w^(-2) + 25/42*w^5 + 25/42*w^(-5);

Besides the obvious w=1, the rational function g-1 has 8 other roots.

> S:= solve(g=1);

These all have absolute value 1.  The corresponding values of x in the interval [0, 10*b] are

> map(t -> 5*b/(I*Pi)*ln(t), {S});

I'm guessing that .9999999996 is really supposed to be 1,  1.380952380 is 29/21, 1.166666667*10^(-14) is 7/(10*b).

Then your expression L[1] is

> f := cos(7*Pi*x/(10*b))*cos(3*Pi*x/(10*b)) - 29/21 * sin(7*Pi*x/(10*b))*sin(3*Pi*x/(10*b));

> combine(f);

-4*cos(2*Pi*x/(5*b))*(1/21)+25*cos(Pi*x/b)*(1/21)

Now if w = exp(I*Pi*x/(5*b)), this can be written as

> g := -2/21*w^2 - 2/21*w^(-2) + 25/42*w^5 + 25/42*w^(-5);

Besides the obvious w=1, the rational function g-1 has 8 other roots.

> S:= solve(g=1);

These all have absolute value 1.  The corresponding values of x in the interval [0, 10*b] are

> map(t -> 5*b/(I*Pi)*ln(t), {S});

Well then, you're asking the wrong question.  You don't need an expression for y(r), you need to know how to calculate integrals involving y(r).  For example, suppose you want to calculate int(y(r)*cos(r), t = 1 .. 2) where y(r) is a solution of differential equation de with initial conditions ic:

> S := dsolve({de, ic}, y(r), numeric, output=listprocedure);
   Y:= subs(S, y(r));

Now the following might work:

> evalf(Int(Y(r)^2 * cos(r), r = 1 .. 2));

If not, you might try this:

> S2:= dsolve({diff(J(r),r) = Y(r)^2 * cos(r), J(1) = 0}, J(r), numeric);
   subs(S2(2), J(r));

Well then, you're asking the wrong question.  You don't need an expression for y(r), you need to know how to calculate integrals involving y(r).  For example, suppose you want to calculate int(y(r)*cos(r), t = 1 .. 2) where y(r) is a solution of differential equation de with initial conditions ic:

> S := dsolve({de, ic}, y(r), numeric, output=listprocedure);
   Y:= subs(S, y(r));

Now the following might work:

> evalf(Int(Y(r)^2 * cos(r), r = 1 .. 2));

If not, you might try this:

> S2:= dsolve({diff(J(r),r) = Y(r)^2 * cos(r), J(1) = 0}, J(r), numeric);
   subs(S2(2), J(r));

OK, now you may have a doable problem.

With the initial conditions y(0)=0, D(y)(0) = 0, even though the DE is singular there, Maple is able to find a series solution.

> de:= diff(y(r),r,r)+2*diff(y(r),r)/r+486*polylog(3/2,-exp(1-1/43*y(r))) = 0;
  ys:= dsolve({de, y(0) = 0, D(y)(0) = 0}, y(r), series, order=10);

ys:=y(r) = (-81*polylog(3/2,-exp(1)))*r^2+(-19683/430*polylog(1/2,-exp(1))*polylog(3/2,-exp(1)))*r^4+(-531441/25886*polylog(-1/2,-exp(1))*polylog(3/2,-exp(1))^2-1594323/129430*polylog(1/2,-exp(1))^2*polylog(3/2,-exp(1)))*r^6+(-186535791/11130980*polylog(-1/2,-exp(1))*polylog(3/2,-exp(1))^2*polylog(1/2,-exp(1))-4782969/636056*polylog(-3/2,-exp(1))*polylog(3/2,-exp(1))^3-43046721/22261960*polylog(1/2,-exp(1))^3*polylog(3/2,-exp(1)))*r^8+O(r^10)

> yn := convert(rhs(evalf(ys)),polynom);

 yn := 127.6269029*r^2-74.07550960*r^4+43.19136809*r^6-22.65830305*r^8

This series appears to have radius of convergence at least 1.  To play it safe, I'll use this for 0 < r < 0.1.  At r = 0.1 we have:

> yn1:= y(0.1) = eval(yn, r = 0.1);
  dyn1:= D(y)(0.1) = eval(diff(yn,r), r = 0.1

 For larger values of r, we may use these as initial values for a numerical solution.

> S:= dsolve({de, yn1, dyn1}, y(r), numeric);

> with(plots): odeplot(S, [r, y(r)], r = 0.01 .. 100);

 For large r, this looks approximately logarithmic.  In fact, I suspect that all solutions are asymptotic to 86*ln(r) + 43 - 43*ln(43/243) as r -> infinity

OK, now you may have a doable problem.

With the initial conditions y(0)=0, D(y)(0) = 0, even though the DE is singular there, Maple is able to find a series solution.

> de:= diff(y(r),r,r)+2*diff(y(r),r)/r+486*polylog(3/2,-exp(1-1/43*y(r))) = 0;
  ys:= dsolve({de, y(0) = 0, D(y)(0) = 0}, y(r), series, order=10);

ys:=y(r) = (-81*polylog(3/2,-exp(1)))*r^2+(-19683/430*polylog(1/2,-exp(1))*polylog(3/2,-exp(1)))*r^4+(-531441/25886*polylog(-1/2,-exp(1))*polylog(3/2,-exp(1))^2-1594323/129430*polylog(1/2,-exp(1))^2*polylog(3/2,-exp(1)))*r^6+(-186535791/11130980*polylog(-1/2,-exp(1))*polylog(3/2,-exp(1))^2*polylog(1/2,-exp(1))-4782969/636056*polylog(-3/2,-exp(1))*polylog(3/2,-exp(1))^3-43046721/22261960*polylog(1/2,-exp(1))^3*polylog(3/2,-exp(1)))*r^8+O(r^10)

> yn := convert(rhs(evalf(ys)),polynom);

 yn := 127.6269029*r^2-74.07550960*r^4+43.19136809*r^6-22.65830305*r^8

This series appears to have radius of convergence at least 1.  To play it safe, I'll use this for 0 < r < 0.1.  At r = 0.1 we have:

> yn1:= y(0.1) = eval(yn, r = 0.1);
  dyn1:= D(y)(0.1) = eval(diff(yn,r), r = 0.1

 For larger values of r, we may use these as initial values for a numerical solution.

> S:= dsolve({de, yn1, dyn1}, y(r), numeric);

> with(plots): odeplot(S, [r, y(r)], r = 0.01 .. 100);

 For large r, this looks approximately logarithmic.  In fact, I suspect that all solutions are asymptotic to 86*ln(r) + 43 - 43*ln(43/243) as r -> infinity

What do you mean by the phaseplane of two variables in a 4-variable system?  The "arrows" would have to depend on the other two variables.

What do you mean by the phaseplane of two variables in a 4-variable system?  The "arrows" would have to depend on the other two variables.

No, it doesn't clarify anything.  I was not missing the negative sign.  The function polylog(3/2, z) is increasing on the interval [-1,1]  and has no real roots other than z=0. 

No, it doesn't clarify anything.  I was not missing the negative sign.  The function polylog(3/2, z) is increasing on the interval [-1,1]  and has no real roots other than z=0.