There are no contradictions. The issue that you are plotting the wrong surface for which you want to obtain contours. Should be

f := exp(-r/2)*r^2*cos(2*phi);
plot3d(f, phi = 0..2*Pi, r = 0..15,  coords = cylindrical);
plots:-contourplot(f,  phi = 0..2*Pi, r = 0..15, coords = cylindrical, grid=[100,100]);
                                

Please note where there is a mark 15.

restart;
sol:=dsolve({diff(ca(t), t) = -3.600000000*10^20*exp(-15098.13790/(340-20*ca(t)))*ca(t), ca(0) = 2}, numeric):
f:=x->eval(ca(t), sol(x));
fsolve('f(t)'=0.2);  # The root
plots:-display(plot(0.2, x=0..5, linestyle=3, color=black),plots:-odeplot(sol,[t,ca(t)], t=0..5, color=red));  # Visual check

restart;
a1 := -k*L*(theta-1)*((z-1)*epsilon*delta-epsilon*z+z);
e:=op([-1,1], a1);
subs(e = freeze(e), a1);
applyop(collect, -1, %, z);
b1:=thaw(%);
                            
 

n:=4:
R:=Matrix(1,n):
for i from 1 to n do R[1,i]:=unapply(i*x, x);
end do:
R;

And what do you expect as a plot? Even if you have 2 parameters  t  and  s , and the functions  x(t,s)  and  y(t,s)  then the set of points on the plane  (x(t,s), y(t,s))  is not a curve, but some flat region. Here is an example of constructing such a region:

plot([seq([t-s,t^2+s, t=0..1], s=0..1, 0.01)], thickness=4, color=green);

I used  plots:-logplot  so that you could view all the curves on one plot:

restart;
a[0], a[1], a[2] := 3.5927*10^33, 1.3197*10^19, 4.8478*10^4: 
A := plots:-logplot([seq(a[i]*x^i, i = 0 .. 2)], x = 2.7*10^14 .. 5.4*10^14, color = [red, blue, green]):
a[0], a[1], a[2] := 2.6235*10^34, 4.876*10^19, 9.0612*10^4:
B := plots:-logplot([seq(a[i]*x^i, i = 0 .. 2)], x = 5.4*10^14 .. 8.2450*10^14, color = [red, blue, green]):
a[0], a[1], a[2] := 8.5561*10^34, 1.0377*10^20, 1.3197*10^19: 
C := plots:-logplot([seq(a[i]*x^i, i = 0 .. 2)], x = 8.2450*10^14 .. 1.80*10^15, color = [red, blue, green]): 
plots:-display(A, B, C);

A:=Array(0..3, 0..3, (i,j)->4*i+j+1);
convert(A, Matrix);

Example for an arbitrary list:

L:=[seq(rand(0..9)(), i=1..16)];
A:=Array(0..3, 0..3, (i,j)->L[4*i+j+1]);
convert(A, Matrix);

f := x->piecewise(1 <= x and x <= 2, c[1]*x+c[2], 2 <= x and x <= 3, c[3]*x+c[4], 3 <= x and x <= 4, c[5]*x+c[6]) ):
op(4, f(x));

L:=[[761, 768, 776, 784, 793, 803, 813, 823, 833, 842], [723, 725, 728, 731, 734, 738, 743, 749, 756, 764], [516, 516, 516, 517, 519, 522, 526, 531, 536, 541], [382, 384, 386, 389, 393, 398, 404, 411, 419, 427], [78, 86, 95, 105, 115, 125, 135, 144, 154, 164], [751, 760, 770, 780, 790, 799, 809, 819, 829, 839], [773, 783, 793, 803, 812, 822, 831, 840, 850, 859], [160, 170, 180, 189, 199, 209, 219, 229, 239, 249]]:
X:=[$ 1..10]:
P:=map(t->`[]`~(X, t), L);
plot(P, color=[red,blue,green,yellow,cyan,gold,pink,violet], size=[750,350], labels=[time,positions]);

We denote  a=-(epsilon-1)*psi/epsilon :

e1:=n = (theta-1)*(z-1)*(-k*psi*L*(-(epsilon-1)*psi/epsilon)^epsilon*k^epsilon+K*(-(epsilon-1)*psi*k/epsilon)^epsilon)/(gamma*((theta-1)*(z-1)*(-(epsilon-1)*psi*k/epsilon)^epsilon-(-(epsilon-1)*psi/epsilon)^epsilon*k^epsilon*theta*z*psi));
    subs(epsilon-1=-a*epsilon/psi, e1);
simplify(%) assuming a>0, k>0;

with(plottools):
A:=rectangle([0, 0], [10, 10], style=line, thickness=2,color = red):
f:=p->rotate(p, arctan(1/3), [5,5]):
g:=p->homothety(p, sqrt(10)/4, [5,5]):
h:=g@f:
plots:-display(seq((h@@n)(A), n=0..10));

It looks like a bug.

3 different results. The plot confirms the correctness of the latter:

z:=exp(I*t):
f:=z^(1/2)*diff(z, t):
int(f, t=0..2*Pi);
int(evalc(f), t=0..2*Pi);
int(simplify(evalc(f)), t=0..2*Pi);
plot([Re(f), Im(f)], t=0..2*Pi, color=[red,blue], scaling=constrained);

Addition. For the second example, all methods (in Maple 2016.2) yield the same result:

f:=1/2+exp(I*t):
evalf(Int(f, t=0..2*Pi));
int(f, t=0..2*Pi);
int(evalc(f), t=0..2*Pi);
int(simplify(evalc(f)), t=0..2*Pi);
plot([Re(f), Im(f)], t=0..2*Pi, color=[red,blue], scaling=constrained);
 

Use  Array  and  Matrix  instead  array  and  matrix:

U := Array(1 .. 5, 1 .. 5);
convert(U, Matrix);
for i to 5 do
for j to 5 do
if i = j then U[i, j] := 1 else
U[i, j] := 0; 
end if end do end do: 
U;

A more compact ways for solution the same task:

U := Array(1 .. 5, 1 .. 5, (i, j) ->`if`(i = j, 1, 0));

or

U := Matrix(5,  (i, j)-> `if`(i = j, 1, 0));

Just assigning  Fh1  and  theta  values from the first solution:

sol := {Fh1 = 121.6477702, theta = 0.9606764638}, {Fh1 = -121.6477702, theta = -2.180916190}:
assign(sol[1]);
Fh1, theta; 
                                         121.6477702, 0.9606764638
 

If you want to save the obtained points  (x,y)  in a list to use later, for example, for plotting the function  y(x), you can assign the resulting values to indexed variables  x[i], y[i]  at each step of the loop:

restart;
x[0]:=0: y[0]:=1: xf:=1: n:=10:
h:=evalf((xf-x[0])/n):
f:=(x,y) -> x+y:
for i from 1 to n do 
    k:=f(x[i-1],y[i-1]);
    y[i]:=y[i-1]+h*k;
    x[i]:=x[i-1]+h;
od:
L:=[seq([x[i],y[i]], i=0..n)];   # List of received points
dsolve({diff(g(t),t)=t+g(t), g(0)=1}, g(t));  # The exact solution
assign(%):
plot([L, g(t)], t=0..1, color=[red, blue],  labels=["x", "y(x)"]);  # The plots of the approximate (red) and the exact (blue) solutions