For Your matrix you can use this:

F := (n, f) -> LinearAlgebra:-BandMatrix(map(f, [$1 .. n]), 0, n);
F(4, (n) -> d[n]);

Thumb if You like.

for i from 1 to n do proc(i, a, b) a[i]:=b[i] end(i, a ,b) od;

also

for i from 1 to n do proc(i) a[i]:=b[i] end(i) od;

provided a, b are not declared local.

Your formula appears to be:

fgP,fgQ:=unapply(gP,x,y), unapply(gQ,x,y);
(fgP(QS1,QS2)*fgQ(S1,-S2)/(fgP(PS1,PS2)*fgQ(S1,-S2))) mod p;

Notice the cancellation.

e_n_1b := n_1 = (-w^(2*sigma)*tau + w^sigma)*s*nsp_1/(w^sigma*tau^2 - w^(2*sigma)*tau + w^sigma - tau) + tau*(w^(sigma - 1)*tau - w^(2*sigma - 1))*s*nsp_2/(w^sigma*tau^2 - w^(2*sigma)*tau + w^sigma - tau);

wpowers := (A, e) -> frontend(convert, [subsindets(A, dependent(w)^anything, expand), parfrac, e], [{Non}(identical(w^sigma)), {}]);
wpowers(e_n_1b, w);
wpowers(e_n_1b, w^sigma);

Thumb if You like.

Edit: There was an error in my response, just like there was one in your code.
Correction:

PlotGraph := proc(func::anything, opts::(seq({identical(zoom) = [range, range]})))
   plots:-display(plot(func), ifelse(membertype(identical(zoom) = [range, range], [opts]), view = subs([opts], zoom), 'NULL'))
end proc;;
PlotGraph(x^2, zoom = [0 .. 1, -5 .. 5]);

Thumb if You like.

One variable, constant bounds:

a := solve({-1 < x, 0 < (2*x)/(x^2 - 1), x < 1}, {x});
solve(`and`(op(a)));

Thumb if You like.

See this example:

A := map((x -> ifelse( x<1/3 or x>2/3, x, 0)), LinearAlgebra:-RandomMatrix(25, 25, generator = rand(0 .. 1.0))):
P := plots[matrixplot](A, heights = histogram, transparency = 0.3, style = patchnogrid, shading = zhue, orientation=[30, 30, 0]):
subsindets(P, [[anything, anything, 0.]$4], (x -> NULL));

Thumb if You like.

The point at K=4*I..1+5*I that appears a solution is at:

{K = (1/2)*sqrt(nu*Pi*(-9*Pi*nu+(8*I)*h^2))/nu};

The above is a solution to cosh(Q)=0. cosh(Q) appears in the denominator.

Thumb if You like.

A := (4*x + 1)/(diff(f(y), y)) = a;
b := (3*(4*x + 1))/((diff(f(y), y))*(3*x + 1));
frontend(algsubs, [A, b, diff(f(y), y)], [{Non}(function), {}]); 

Thumb if You like.

Do You mean to do something like this:
"numer(L = 2/3)";
f := e -> lcm((denom@lhs, denom@rhs)(e))*e;
f(L = 2/3);

 

collect(Your_stuff_goes_here, s, factor); #normal might do also

Thumb if You like.

You can plot this directly by applying the transform to the original inequalities:

  T[0] := [abs(z) < 3, 1 < abs(z - 1)];
  T[1] := z = x + I*y;
  T[2] := evalc(subs(T[1], T[0])) assuming real;
  T[3] := w = radnormal(subs(x1 = 9/2 - 1/2*sqrt(45), x2 = 9/2 + 1/2*sqrt(45), ((3 + sqrt(5))*(z - x1))/(2*(z - x2))));
  T[4] := z = solve(T[3], z);
  T[5] := radnormal(subs(T[4], T[0]));
  T[6] := w = I*v + u;
  T[7] := factor(expand(subs(T[6], T[5]))) assuming real;
  plots[inequal](T[2], x = -4 .. 4, y = -4 .. 4);
  plots[inequal](T[7], u = -3/2 .. 3/2, v = -3/2 .. 3/2);

Thumb if You like.

Look at:

Sol := [solve(x^3+(a-3)^3*x^2-x*a^2+a^3 = 0, x, real, parametric)];

Also consider this method:

Sol2 := [solve(x^3 + (a - 3)^3*x^2 - a^2*x + a^3 = 0, a, real, parametric)];
subs(Sol2[1][1], a);
plot(%, x = -4 .. 4);
#or
#Sol2 := [RealDomain:-solve(x^3 + (a - 3)^3*x^2 - a^2*x + a^3 = 0, {a})];
#subs(Sol2[1], a);
#plot(%[1], x = -4 .. 4);

Thumb if You like.

Notice your E procedure is missing x as input.

Simply use:
(f@@depth)(x);

You have to define the region of integration in a new way:

Int(Int(Int(1, y = 0 .. 1), tau = 0 .. z + q), z = 0 .. t); 
IntegrationTools[CollapseNested](%); 
subs([(tau = 0 .. z + q) = (z = piecewise(tau < q, 0, q <= tau, tau - q) .. t), (z = 0 .. t) = (tau = 0 .. q + t)], %); 
(expand(value(%%) = value(%)) assuming (0 < t, 0 < q)), (value(%%) = value(%));

Thumb if You like.