I think that this is the case when it is easier to solve by hand, because it is easy to express exp(x) rationally through tanh(x/2). If you need to perform a similar conversion more than once, the procedure called  Convert  will be useful. It converts any occurrence of  exp  function into  tanh  function:

restart;
Convert:=proc(Expr)
local A;
convert(tanh(z), exp);
combine(normal(subs(exp(-z)=1/exp(z), %)));
solve(%=tanh(z), exp(2*z));
A:=subs(z=z/2, %);
simplify(applyrule(exp(z::anything)=A, Expr));
end proc:

Examples of use:

Convert(exp(x));
Expr:=-(exp(-a*s)-1)*kw/((exp(-a*s)+1)*s^2);
Convert(Expr);  # The initial example

Maple uses radians by default, so should be Pi/2 instead of 90.

Output := Array(-10 .. 10, -10 .. 10, (i, j)->i^2+j^2):
M:=convert(Output, listlist):
plots:-listcontplot(M, axes=none);

lambda:=2.6:
fsolve(sum(lambda^k*exp(-lambda)/k!, k=0..n)=0.95);

                                  4.980237188

Example:

restart;
sol:=solve({x-2*y=5, x+3*y=1});
eval(x, sol);  # The first way
eval(y, sol);
assign(sol);  # The second way
x, y;

If you need decimals, use  fsolve  instead of  solve.

Temporarily replace symbols by unique numeric constants.

Example:
P:=a*x^5+a*y^5-b*x*y+c:
Q:=subs([a=Pi, b=exp(1), c=gamma], P);
convert(Q, elsymfun);
subs([Pi=a, exp(1)=b, gamma=c], %);

I did not understand how you define matrix  A  for your example, I used  rand  command for this.

restart;
lavplot := proc() 
local i, j, feltfarver, points, colors, x, y; 
global A;
uses LinearAlgebra, plots; 
x := ColumnDimension(A); 
y := RowDimension(A); 
points := []; 
feltfarver := [blue, green, "PeachPuff", "DarkGreen", "Gold", "Silver", "Chocolate", yellow]; 
colors := []; 
for i to y do 
for j to `if`(i::odd, x, x-1) do 
if A[i, j] <> 0 then points := [op(points), [`if`(i::odd, j, j+1/2), (x+1-i)*sqrt(3)/2]]; colors := [op(colors), feltfarver[A[i, j]]] end if; 
end do; end do; 
pointplot(points, symbol = solidcircle, color = colors, symbolsize = 60, axes = none, scaling=constrained); 
end proc:

Example of use:

r := rand(1 .. 8): 
A := Matrix(7, 10, (i, j)->r()); 
lavplot();

Edit.

The same problem can easily be solved with the help of  plots:-animate  command. This is more convenient, because it allows you to easily adjust the number of frames and the step of the animation. In the animation panel, you can see the same slider, and you can also sequentially output frames one by one in any direction. For example, below we choose the range for the parameter a = -2..2  and animation's step=0.1 and then we have to take the number of frames  (2-(-2))/0.1+1=41

restart:
with(plots):
animate(inequal,[{x>0,x>2*(y-a-1),y>x+a-1}, x=-1..4, y=-3..5, scaling=constrained], a=-2..2, frames=41);

Clearly, the solution is the interior of a triangle as the intersection of three half-planes. In the system, I removed the inequality  x+a+1>y, because for x>0 we have  x>2*(y-a-1) => x+a+1>y

Edit.

I do not know of a simple method that would lead to the result as in a textbook (or a close result). Here is a way (partially manual), which in 3 steps gives the result almost as in a textbook:

A:=int(x*sqrt(a*x^2+b*x+c), x);

B:=simplify(A, symbolic); # First step
C:=`*`(op(1..2,B))*map(normal,`+`(op([3,1..2],B))); # Second step
`*`(op(1..2,C))/(-12)*``(op(3,C)*(-12));  # Result

The result:

For details on numerical integration see Help  ?evalf/int

 idol050279 10
Inside the opaque cylinder, 30 random points are located. When the plane moves, the walls of the cylinder to the left of the plane disappear and the points become visible:

restart;
with(plots): with(plottools):
P1:=[seq(RandomTools:-Generate(list(float(range = -3..3), 3)), n=1..300)]:
P2:=select(p->p[1]^2+p[2]^2<3.9 and p[3]>0, P1)[1..30]:
Al:=animate(plot3d,[[2*cos(t),2*sin(t),h], t=arccos((2-s)/2)..2*Pi-arccos((2-s)/2), h=0..3, style=surface], s=0..4, frames=90, paraminfo=false):
Ac:=animate(plot3d,[[[x,y,3],[x,y,0]],x=2-s..-2,y=-sqrt(4-x^2)..sqrt(4-x^2), style=surface], s=0..4, frames=90, paraminfo=false):
B:=plots:-pointplot3d(P2, symbol=solidsphere, color=red, symbolsize=10):
C:=animate(display@polygon,[[[2-s,2,0],[2-s,2,3],[2-s,-2,3],[2-s,-2,0]], color=yellow], s=-0.02..4, frames=90, paraminfo=false):
plots:-display(Al, Ac, B, C, axes=normal, view=[-2.7..2.7,-2.7..2.7,-0.7..3.7]);

Example:

msolve(2^n=2, 3);
                                            {n = 1 + 2 _Z2}

_Z2  is an integer.
 

You forgot to put a semicolon in 1dmath mode. Also before  @  you have confused the order  "  and  ':

["<","&",">","'@"];

                                            ["<","&",">",'"@"]

Alternative solution:

restart;
f:=1-max(abs(x),abs(y),abs(z)):
g:=1.5-(x^2+y^2+z^2):
A:=plots:-implicitplot3d(piecewise(f<g,f,g), x=-1.6..1.6,y=-1.6..1.6, z=-1.6..1.6, 
       numpoints=300000, scaling=constrained, color=green, style=surface, axes=none ):
f:=1-(x^2+y^2+z^2):
g:=z^2-x^2-y^2:
B:=plots:-implicitplot3d(piecewise(f<g,f,g), x=-1..1,y=-1..1, z=-1..1, 
       numpoints=300000, scaling=constrained, color=green, style=surface, axes=none):
plots:-display(<A | B>);

Interestingly, for a slightly smaller numpoints, the plotting quality is slightly better.

See this simple example in which the difference between  algsubs  and  subs  is clearly visible. Your second example works as expected, because the numerator is independent of  w :

Expr:=a/(a+b);
algsubs(a+b=c, Expr);
subs(a+b=c, Expr);

                         Expr := a/(a+b)
                               (-b+c)/c
                                  a/c

Read Help on these commands for details.