restart;
myPi_1:=proc(r)
local i,l;
l:=0;
for i from 1 to floor(r) do
    if isprime(i) and (i mod 4 = 1) then
        l:=l+1;
    end if;
end do;
return l;
end proc:    

plot('myPi_1(t)', t = 0 .. 500, numpoints = 500, thickness = 2, color = black);

To calculate the values of functions  U1  and  U, it is better to define these functions as procedures and calculate their values for each call of these procedures. Then you will not have any problems:

G1:=-0.9445379894:
f:= (x) -> 0.9/abs(x-0.4)^(1/3)+0.1/abs(x-0.6)^(1/2):
U1 := x->-exp(-x)*(int(f(t)*exp(t), t = 0 .. x)+G1)/2-exp(x)*(int(f(t)*exp(-t), t = 0 .. x)+G1)/2;
seq(U1(x), x=0..1.5, 0.1);  # Calculations of values the function U1(x)  in the different points
plot(U1, 0..1);  # The plot of U1(x)

For a function  U , everything is the same.

lst:=[1,2,3,4,5,6]:
lst[2..5];
                                                        [2, 3, 4, 5]

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