Use the command odetest.

First note that in differential equations, derivatives must be entered as diff(E(r),r), not diff(E,r). Here's your system of ODEs, corrected:

odesys:= {
   diff(E(r),r)^2+diff(B(r),r)/E(r)+E(r)*B(r)*r=0,
   diff(B(r),r)^2/E(r)+5*r = 0
}:

And here's your putative solution:

puta:= {E(r) = r^a, B(r) = exp(a*r)}:

odetest(puta, odesys);

The fact that the above is not identically (with respect to r) zero for any a shows that puta is not a true solution.

I only have time right now to give a brief example of how this might be accomplished.

The interpolation will be handled by command CurveFitting:-ArrayInterpolation (see ?CurveFitting,ArrayInterpolation). In the example, I construct arrays that will interpolate an approximation to cosine. Then I write a procedure that does the call to ArrayInterpolation. Then I use this procedure in a differential equation whose solution should be an approximation of sine. The interpolating procedure must be declared to dsolve as known.

Download known.mw

@Markiyan Hirnyk 

To avoid catching expressions with integer exponents, you need a more-specific type specification. This works for any expression ex:

indets(ex, algebraic^fraction)

The sqrt(-1), aka I, needs to be treated separately. If you want to catch this also, use

indets(ex, {identical(sqrt(-1)), algebraic^fraction}),

which will work even if the imaginary unit has been redefined.

Complicated roots are sometimes expressed in RootOf form. To catch these also, use

indets(ex, {specfunc(polynom, RootOf), identical(sqrt(-1)), algebraic^fraction}).

Here's something to try; I don't know if it'll work. Open the file in Maple 16, and then save the expression using the command

save expr "s.maple.m";

The .m causes the file to be saved in Maple internal format. This may be easier for Maple 18 to read.

What I see from a quick read through is that you never define r1 or r2. You can't plot them if you don't define them.

Here's a basic module factory to accumulate arithmetic means (not exponential means). I just wanted to show the basics of a module factory, and how simple they can be.

Download module_factory.mw

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"])