@MuriloLobo Something is missing. As Mac Dude said, S1 is not defined anywhere in there. You should be getting Error, `S1` does not evaluate to a module. You also need a semicolon at the end of each line. The comma at the end of the last seq needs to be a semicolon, as Mac Dude said.

Please upload your worksheet, don't just post it. Use the green uparrow that is the last icon of the second row of the toolbar of trhe MaplePrimes editor.

@Mac Dude Creating a sequence of equations or inequalities with seq is fine. Nearly every structure in Maple (everything that I can think of) is a first-class object, i.e., it can be assigned to a variable and manipulated like any basic variable.

@radzys Everything works in the worksheet that you posted: The solve solves the eight equations and you get the animation at the end. If you still have an issue that you'd like me to look at (your "most perfect scenario"), please post a worksheet where you attempt it and it doesn't work.

The period of the oscillations is about 0.0006. This is too small for you to see any significant detail on the scale x=0..1. Better is x=0..0.001.

@OV I am referring to calling s1[3](1) after running the loop as modified by Mehdi and getting the recursion error in "_dat".

There is nothing wrong with having an Array of procedures in Maple. (Technically, this is a table of procedures, but the distinction between table and Array is insiginificant here.) The problem is that s1[3], which is a procedure generated by dsolve, calls itself. I still don't know how that happened.

@Michael eval instead of evalf will produce exact results. But better still to use two-argument eval, as Acer did.

@Alejandro Jakubi The controller gets swapped out whenever it is in a waiting state, which is usually most of the time. This leaves all cpus available to run grid kernels.

Please post your worksheet.

@mehdi jafari The second recursion error is in what the OP calls s1[3], which is a procedure extracted from dsolve(..., numeric) output. It is in the second cycle of the for that this error occurs. It is confusing, and perhaps the source of the error, that s1[...] is defined at two places within the loop with different indices.

You can cause the error simply by calling s1[3](1) or at any numeric input.

@Markiyan Hirnyk Yes, all the necessary assuming is done by evalc in that case. I made a mistake by noticing that evalc does automatically map over lists, and then failing to notice that it does not also automatically map over vectors.

@rashmi Show me exactly what you did to get plunging orbits.

@Alejandro Jakubi I don't like that my method requires two sets of quotes; I can't explain that. But note that the % method requires the use of % in the target expression as well as the rule. This means that a substitution needs to be made in the OP's Intat expression since it was programmatically generated. Also, the % technique will ultimately require two applications of value: the first will change %Intat to Intat (and change Sum to sum), and the second will change Intat to intat. To overcome this, a procedure can be made

ToSumIntat:= (J::'Intat'(algebraic, symbol= algebraic))->
     value(
          applyrule(
               %Intat(
                    A::algebraic*Sum(v::algebraic,r::{name,equation})*
                    B::algebraic,u::equation
               )= %Sum(%Intat(A*v*B,u),r),
               subs(Intat= %Intat, J)
          )
     )
;

Note the use of %Sum on the rhs of the rule.

So, I think that the % technique is ultimately more to type.

Please give us the PDE and the pdsolve command that you used to get the Whittaker functions.

Give a specific example please.

@rashmi Any plots on the same axes can be combined with plots:-display. Or you could specify multiple plots in the odeplot command, like this:

Sol:= dsolve(
     {diff(r1(t),t)=cos(t), r1(0)=0, r2(t)=2*sin(t), r2(0)=0},
     {r1(t), r2(t)},
     numeric
):
plots:-changecoords(
     plots:-odeplot(Sol, [[r1(t), t], [r2(t),t]], axiscoordinates= polar),
     polar
);