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

You have an equation of the 9th degree with respect to lambda with two parameters  x  and  y . It is well known that in the general case the roots of such an equation can not be expressed in principle in terms of its coefficients. You can solve this equation numerically, but only by first setting the parameters.

Example:

x:=1: y:=2:
P:=_Z^9-28*_Z^8+328*_Z^7-2088*_Z^6+(-2*cos(y)-2*cos(x)+7852)*_Z^5+(28*cos(y)+28*cos(x)-17768)*_Z^4+(-152*cos(y)-152*cos(x)+23632)*_Z^3+(408*cos(y)+408*cos(x)-17232)*_Z^2+(-12*cos(x)*cos(y)-540*cos(y)-540*cos(x)+5844)*_Z+32*cos(x)*cos(y)+256*cos(y)+256*cos(x)-544:
fsolve(P, _Z);

   0.1348905839, 0.8646312171, 1.143198734, 1.706737179, 3.176390681, 4.304915397, 4.893289594, 5.142364314, 6.633582298

f:= x-> diff(g(x),x)/(1+g(x)):
applyop(eval, 2, f(x), g(x)=h);
subs(h=g(x), series(%, h=0, 4));

I do not have Maple 18 and so I can not confirm it, but try these options:

plot(x^2, x=-2.5..2.5, tickmarks=[spacing(0.5), default]);
# Or   
plot(x^2, x=-2.5..2.5, tickmarks=[[seq(-2.5+0.5*i=-2.5+0.5*i, i=0..10)], default]);