Use  solve  command, which symbolically solves systems of polynomial equations in the form RootOf of a high degree polynomial, after that  allvalues  command and then apply  evalf  command.

Edit.

In your integral,  y  variable is the integration variable (dummy variable). A definite integral does not depend on the name of the variable of integration, so simply replace y under the integral sign with any symbol that is not in your worksheet. Here is a simple example in which the same error is encountered:

f := x->int(x/(x+sin(x)), x = 0 .. x^2); 
f(1/x);   # An error
F := x->int(s/(s+sin(s)), s = 0 .. x^2);  # Workaround
F(1/x);  # OK

Edit.

Example:

randpoly([x, y], dense, degree = 4);
select(t->degree(t,x)+degree(t,y)<2, %);  

Or even shorter:

select(t->degree(t, [x,y])<2, %);

Edit.

Simple procedure Cramer returns the main determinant of the system  det[A] , all the auxiliary determinants  det[i]  and the values of all unknowns calculated according to Cramer's rule:

Cramer:=proc(Sys::{set(equation), list(equation)}, Var::list(name))
local n, A, V, det, S;
uses LinearAlgebra;
n:=nops(Sys);
A, V := GenerateMatrix(Sys, Var);
det[A]:=Determinant(A);
S:=seq(Determinant(<A[..,1..i-1] | V | A[..,i+1..n]>), i=1..n);
print(det__A=det[A], seq(det[i]=S[i], i=1..n));
seq(Var[i]=S[i]/det[A], i=1..n);
end proc:

Example of use:

A:=LinearAlgebra:-RandomMatrix(4, generator=-9..9):
V:=LinearAlgebra:-RandomVector(4, generator=-9..9):
Sys:=LinearAlgebra:-GenerateEquations(A, [seq(x[i], i=1..4)], V);
Cramer(Sys, [seq(x[i], i=1..4)]);
    
 

restart;
with(DEtools): 
DE1 := diff(y(x),x) = diff(y(x),x$3)+3*diff(y(x),x$2)+4*diff(y(x),x)+12*y(x);
DEplot(DE1, y(x), x = 0..5, [[y(0)=-3, D(y)(0)=0, (D@@2)(y)(0)=0]]);

Your equation only Maple 2017 will be able to solve:

pde := (diff(u(x, y), x, x))+diff(u(x, y), y, y) = 0;
 bc[1]:= D[2](u)(x,0) = 0;
 bc[2]:= u(x,1) = x^2-x;
 bc[3]:= u(0,y)=0;
 bc[4]:= u(1,y)=0;
 pdsolve({pde, bc[1], bc[2],bc[3],bc[4]}, HINT= X(x)*Y(y));

HINT  option can be omitted.

X:=<1,2,3,4>:
Y:=<5,6,7,8>:
x.X+y.Y;

In this example  X  and  Y  are vectors. For matrices everything is the same.

Edit.

Try  ContoursWithLabels  procedure  from  here .

Your example:

ContoursWithLabels(-3.392318438*exp(-4.264628895*x)*sin(6.285714286*y), -1/2 .. 1/2, -1/2 .. 3/2, {seq(-12 .. 12, 3)}, [color = black, axes = box, size=[1000,600]], Coloring = [colorstyle = HUE, colorscheme = ["Cyan", "Red"], style = surface]);

If you do not want coloring, then just remove  Coloring  option from the code.

With the built-in command  plots:-contourplot , do so

plots:-contourplot(-3.392318438*exp(-4.264628895*x)*sin(6.285714286*y), x = -1/2 .. 1/2, y = -1/2 .. 3/2, contours=[seq(-12..12,3)], size=[800,500], numpoints=10000);

Edit.

Execute:

op~(Sol);

Do := proc( F::list(procedure) , Q::list(list))
local n:=nops(F), m:=nops(Q);
if n<>m then error "Should be nops(F)=nops(Q)" fi;
seq(F[i]~(Q[i]), i=1..n);
end proc:

Example of use:

Do([x -> x , y -> y^2], [[0,1,2,3], [4,5,6,7]]);

                                            [0, 1, 2, 3], [16, 25, 36, 49]

Edit.

I took  Gamma/(2*Pi)=1  and plotted  Re(w)  (in red) and  Im(w)  (in blue):

restart;
w:=theta-I*ln(r):
X, Y:=(Re,Im)(w) assuming theta::real, r>0;
plot3d([[r*cos(theta),r*sin(theta),X], [r*cos(theta),r*sin(theta),Y]], theta=0..2*Pi, r=0..3, color=[red,blue], axes=normal, view=[-4.3..4.3, -4.3..4.3, -1..6.7]); 

See help on  ?plot,tickmarks

Example:
restart;
Digits:=20:
plot(x^2, x=0..2, 1..1.00000000001);  # Long values
plot(x^2, x=0..2, 1..1.00000000001, tickmarks=[default, [seq(1+2*10^(-12)*k=1+2*k*10^(`-12`), k=1..5)]]);  # Short values

Use nested  seq  command.

Example:

seq(seq(A[i,j] >= 0, j=1..4), i=1..3);

You can not find a limit with your procedure, because it for any particular n simply returns a number, but an exact dependence of the sum on n is needed.
The limit can be found as follows:

S:=unapply(sum(2*i/n*2/n, i=1..n), n);
simplify(S(n));
limit(%, n=infinity);

By  p  procedure the limit can be found only numerically with  a specific accuracy:

evalf[5](p(100000));

                                           2.0000

gl_inveq:=(t,x,alp)->tan((x-t)/alp/(1+t^2))-t:
gl_inv:=(v,alp)->fsolve(gl_inveq(t,v,alp)=0, t):
gl_inv(2,1); # The calculation of a value of the function
evalf(Int(gl_inv(2,t), t=0..1));  # The calculation of the value of the integral

                                            0.8299950485
                                            -0.5405524297