I have a fix and will email it to the poster.  There was an invalid character. I'll update my procedure to handle it.

It works if you use tty Maple (cmaple) rather than Standard Maple to generate the plots.  

You could linearize the equations at a point and then compute the laplace transform of that.  See the help page for DynamicSystems,Linearize.

an := 1/2/Pi*int(x^2*cos(n*x),x=-Pi..Pi) assuming n :: integer;
bn := 1/2/Pi*int(x^2*sin(n*x),x=-Pi..Pi) assuming n :: integer;
a0 := 1/2/Pi*int(x^2*cos(0*x),x=-Pi..Pi);
a0 + sum(an,n=1..infinity);

The following works, but I'm questioning its efficiency. I don't have any better ideas at the moment.

list_minus := proc(L,E)
local R;
    R := Array(E);
    remove(proc(e)
           local p;
               if member(e,R,'p') then
                   R[p] := NULL;
                   true;
               else
                   false;
               end if;
           end proc, L);
end proc:

list_minus([1,2,4,6,2,1,3,6,2],[7,4,2,5,2]) = [1,6,1,3,6,2];

A more compact, but essentially equivalent, implementation is

list_minus := proc(L,E)
local R,p;
    R := Array(E);
    remove( e -> member(e,R,'p') and assign('R[p]') = NULL, L);
end proc:

It would be helpful if you posted a model of interest.  One approach is to use the Linearization template to create a linearized model, then use the LinearAnalysis template to inspect its poles.  

Is this for documentation?  Or do you expect to evaluate the expression in Maple? For documentation you can, in a text region, type C-r to enter inert math, type sum, then use completion (C-shift-space on Linux, not sure otherwise), select the summation symbol with over/under scripts, enter k >= 0 and k <> p in the underscript and the summand in its place.

I don't have a complete, coherent answer.  My surprise is that the call before the restart works. You probably won't be able to use ModuleApply in the fashion you desire. As you undoubtedly know, you can do the following.  Add the static export

Apply :: static := proc( self :: ModularArithmetic, n :: posint )
    modp(n, self:-p);
end proc:

Then use it as

m7 := NewModulo(7):
Apply(m7, 42);
                            0

That will work before and after the restart.

Assuming a is nonnegative, then you want to use the floor.

Carl's suggestion of using an Array of records is a good one, given what we know.  Depending on what you want to do, there could be more efficient (however defined) alternatives. If all the cells have the same set of characteristics, you could assign the finite characteristics to slots in another dimension of the Array. An advantage of that is that you can then cheaply access arbitrary slices of the data.

A record is another possibility.  If you are creating a package with multiple procedures, using records as output is convenient in that the format is partially self-documenting, provided you pick appropriate field names for the record. To return a record you can do

return Record('f' = fval, 'g'=gval);

where fval and gval have been assigned values.

Insert the statement 

k;

after the loop.  If this is part of a larger procedure, you would, instead, use 

print(k);

Better might be to add a message.  One way to do that is to use printf:

printf("value of k is %a\n", k);

Do you want to print them as ascii text, with no 2D math?  In that case you could do something like

val := [0,0,1]:

printf("%d%d%d\n", op(val));
001

If you would rather print the result using 2D output, you could do

nprintf("#mn(\"%d%d%d\");", op(val));
                                         001

You would combine that with the rest of the equation: P(...) = ... ; where the first ellipsis represents the nprintf expression, the second whatever you intend on the rhs.

I used setattribute to save the index of a value with the value, took the minimum of the chosen values, then used attributes to return its attribute, which is the index.

argmin := proc(d :: list, U :: set(posint))
local u;
    attributes(min(seq(setattribute(evalf(d[u]),u), u=U)));
end proc:

argmin([10,7,4,2,5,7],{2,3,6});
                                                  3

You can use the Iterator package in the application center:

with(Iterator):
C := CartesianProduct( L1, L2, L3, ... ):  # L1, L2, etc, are the lists to iterate through
add(f(c), c in C);

If you actually need to iterate through all the values, you can do

for c in C do f(c) end do:

Note that c is assigned a 1D dimensional Array, so f would have to be written accordingly.