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
);

@mehran rajabi 

Setting maxfun= 0 will remove the maxfun restriction.

@Preben Alsholm 

Maple's multiple assignment allows a swap of variables without need of a temporary variable, and it's more efficient to do it that way.

BubbleSort:= proc(L::{list,Array,Vector})
local i, j, T:= Array(L), n:= numelems(L);
     for i to n-1 do
          for j to n-i do
               if T[j] > T[j+1] then
                    T[j],T[j+1]:= T[j+1],T[j] #swap
               end if
          end do
     end do;
     convert(T, whattype(L))
end proc:

BubbleSort(RandomTools:-Generate(list(float, 12)));

@Edinburgh The mistake is your usage of z in two different ways. You have z as the name of the third coordinate function and also as the table that stores the coordinate value. You should make the third coordinate function h to match f and g.