Answering your first question:

local Pi;
cos(x+Pi);

Regarding your second question: I don't know. The answer is not well defined. I mean, why should it return cos(x+Pi) rather then cos(x-Pi)? They're both equal to -cos(x).

There is no need for Arrays as long as you never assign to indexed list entries, which is not allowed for lists longer than 100 elements. By using seq to create lists and map to operate on them, you never need to assign to list entries. Here is a vastly simplified version of your code that produces exactly the same output and works for long lists.

Download plotScript.mw

Your integral equation can be converted to an ODE IVP. Then Maple's numerical ODE solver can be used on it.

Download Int_eqn_to_ode.mw

I am only debugging here the expression which updates z2 because the others were commented out. The errors in the commented-out expressions are very similar.

There are two bugs, both of them in the last factor in your expression for updating z2. First, you have arro[k] immediately in front of the parentheses that surround the (cosine + I*sine) expression. When a name is immediately followed by a left parenthesis, it is interpretted as function application rather than multiplication. You should put an explicit multiplication sign between the arro[k] and the left parenthesis. 

The second bug is that all function applications do require parentheses. That includes sin and cos, even though the parentheses are often omitted in standard mathematical notation and typesetting. So correcting both errors, the final expression may be entered as

arro[k]*(cos(2*Pi*j/21)-I*sin(2*Pi*j/21))

See the shell escape command ?escape.

Here's is a corrected version of your worksheet. I only corrected the errors necessary to get the combined plot. I did not deal with the singularities.

odeplot_(1).mw

The semicolons that you used in the display command should be commas:

 display(p(-1), p(-0.5), p(-0.1));

It's as simple as that.

There're only two uses of the semicolon in Maple: separating complete statements and separating the rows of a Matrix in the abbreviated angle-bracket input.

To create the side-by-side plots, you pass a row vector of plots to plots:-display. The row vector can be created by iterating (with seq) over the values of R exactly like you iterated over the values of J.

eta:= 1+k*x+epsilon*sin(2*Pi*x):
A1:= x-> -exp(-alpha*x)*J^2*R/(2*eta^3+6*xi*J^2*eta^2):
psi0:= A1(x)*y^3:
psi:= delta*psi0:
V:= -diff(psi,x):
delta:= 0.1:
epsilon:= 0.01:
alpha:= 1:
xi:= 0.001:
k:= 0.1:
As:= [0, 2, 4]:
Rs:= [2.0, 5.0, 6.5]:
x:= 0.2:

plots:-display(
     `<|>`(seq(
          plot([seq](
               [eval(V, [J= A, R= r]), y, y= 0..eval(eta, J= A)], A= As),
               title= sprintf("velocity at R=%3.1f", r), labels= ["v", "y"],
               color= [green, red, blue],
               linestyle= [solid, dash, dot],
               legend= [seq](J = A, A= As),
               axes=boxed
          ), r= Rs
     ))
);

The three plots will superficially look alike, but note that the scaling on the horizontal axis is different. I don't know how to post side-by-side plots to MaplePrimes, so just execute the above code and you'll see what I mean.

Using the array form of points[1], give the command

surfdata(convert(points[1], list), dimension= 2, colorscheme= ["Orange","MediumBlue"])

@RAfossey You have invstrans, not invztrans. And why have you separated the two parts of the z expression with a comma? And why did you put y in the command? Where are you getting these wrong ideas from? The command that you need is

invztrans(z/(z-2)*5*z/(5-z)^2, z, n);

That's all that there is to it (for part a)!

For part b) use sum command. Look up ?sum. There's no need for the ztrans command.

Take the answers from parts a and b and form their quotient (as your professor suggests) and apply simplify to it. You should get 1, which shows that the expressions are equal. A more normal way to prove equality is to take the difference of expressions and try to simplify it to 0.

What you want to do can be done simply with piecewise, with no need for events. If the first derivative is 10 or greater, then the second derivative is 0. What did you try with piecewise?

Download piecewise_dsolve.mw

The second parameter of My is y. Inside My you have multiple subexpressions such as diff(a(y), y). This means the derivative of a(y) with respect to (wrt) y. For this to make sense, y has to be a name. But in your call to My, you pass -1/2*b as the second argument. This becomes the value of y. Perhaps what you intended was not the derivative of a(y) wrt y but rather the derivative of a wrt its parameter, then evaluated at y. This would be denoted D(a)(y) rather than diff(a(y),y).

Another problem that I noticed is your use of Pi. You use it both as a constant (as in sin(m*Pi*x/a)) and as a function (as in Pi(y)). It can't be both. Preferably, you should only use it as the constant (its default meaning) and choose another symbol for the Pi in Pi(y).

A similar problem is your usage of m: It is used both as a bound variable (as in m= 1..infinity) and as a function (as in m(y))

Lastly, How can you understand or manage such a long and sloppy function? You should break it up into smaller, more-understandable chunks.

You said that "at best" you'd like to access the points by name. The following lets you do that.

Download Points.mw

The commands matrix, vector, transpose, linsolve, and `&*` are all very old deprecated commands, which means that they shouldn't be used anymore. They only exist in current Maple so that legacy code still works. The `&*` was used for matrix multiplication, and linsolve was used to solve a linear system. In modern Maple the above code is

a:= < a1,a2; b1,b2 >;
b:= < c1,c2 >;
c:= a^%T; #transpose
LinearAlgebra:-LinearSolve(c.a, c.b)

For example:

A_vals:= [1,2,3]:
m:= min(A_vals):  R:= max(A_vals) - m:
Sol:= dsolve(
     {diff(y(x),x$3)=y(x), y(0)=0, D(y)(0)=-0.5, (D@@2)(y)(0)=A},
     parameters= [A], numeric, range= -2..2
):
for v in A_vals do
     Sol(parameters= [v]);
     P[v]:= plots:-odeplot(
          Sol,
          [x,diff(y(x),x)],
          color= COLOR(HUE, evalf(.85*(v-m)/R)),
          legend= [A=v]
     )
end do:
plots:-display(convert(P,set));