I usually use the following way:

restart;
u:=<-2+m,3+m>;

u:=eval(u, m=5);
# or
u1:=eval(u, m=5);

Such equations are easily reducible to ordinary differential equations with derivatives. In the example below, we solve the equation  p*dx-q*dy=1/2*d(y^2-x^2)  by taking  p=3, q=2  and plotting several integral curves and 2 asymptotes:

restart;
d:=z->diff(z,x)*dx:
p:=3: q:=2:
p*d(y(x))-q*dx=1/2*d(y(x)^2-x^2);
dsolve(%);
plot([seq(eval(y(x),%[1]),_C1=-8..3),seq(eval(y(x),%[2]),_C1=-8..3),-(x-2)+3,x-2+3], x=-4..8, color=[green$24,red$2], scaling=constrained);

restart;
S := proc(n) 
local i, s;
uses NumberTheory;
s:=0; 
for i to n do 
if i in Divisors(n) and type(n/i, odd) then s:=s+i end if; 
end do;
s;
end proc:

# Example
S(12);

                                                  16

The code below can be made more automatic, but it's a rather tedious job. Let it be for now:

restart;
expr1:=sqrt(x^2*y - 4)*sqrt(x^2*y);
expr2:=sqrt((x^2*y - 4)*(x^2*y));
SolveTools:-SemiAlgebraic({op([1,1],expr1)>=0,op([2,1],expr1)>=0, op([1,1],expr2)>=0 , expr1^2=expr2^2}, {x,y});

restart;

plot3d(y^2, x=0..1, y=-1..1, style=surface, color=grey, filled=true, axes=normal, view=[-0.5..1.5,-1.5..1.5,0..1.5]);

Volume=int(y^2, [y=-1..1,x=0..1]);

plots:-implicitplot(x=y^2, x=0..9, y=-3..3);

The programm uses the first  2*K terms of your list only (you have  K=1 ). So all works properly:

restart;
interface(version);

genfunc:-rgf_findrecur(1, [1, 2], t, n);
t:=unapply(rsolve({%, t(1)=1},t(n)), n);
seq(t(n), n=1..4);

The best quality will be if all surfaces are parameterized:

restart;
x:=r*cos(phi); y:=r*sin(phi); z:=3-x; H:=eval(z,r=2);
S1:=plot3d([x,y,0], r=0..2, phi=0..2*Pi, style=surface, color=blue): # Bottom
S2:=plot3d([eval([x,y],r=2)[],h], phi=0..2*Pi, h=0..H, style=surface, color=green): # Side surface
S3:=plot3d([x,y,z], r=0..2, phi=0..2*Pi, style=surface, color=red): # Top surface
plots:-display(S1,S2,S3, axes=normal, view=[-3..3,-3..3,0..5], scaling=constrained, lightmodel=light4, labels=["x","y","z"], labelfont=[times,bold,14]);

We can also show the cutting plane using the  transparency  option:

restart;
x:=r*cos(phi); y:=r*sin(phi); z:=3-x; H:=eval(z,r=2);
S1:=plot3d([x,y,0], r=0..2, phi=0..2*Pi, style=surface, color=blue):
S2:=plot3d([eval([x,y],r=2)[],h], phi=0..2*Pi, h=0..H, style=surface, color=green):
S3:=plot3d([x,y,z], r=0..2, phi=0..2*Pi,  style=surface, color=red):
S4:=plot3d(3-u, u=-3..3, v=-3..3, style=surface, color=yellow, transparency=0.5):
plots:-display(S1,S2,S3,S4, axes=normal, view=[-3..3,-3..3,0..5], scaling=constrained, lightmodel=light4, labels=["x","y","z"], labelfont=[times,bold,14]);

                          

Edit.

f:=1/(z^2+I)^2;
convert(f, parfrac, sqrt(2));

We obtain this number  sqrt(2)  (a field expansion) by solving the equation

solve(z^2+I, z);

Some expressions are written ambiguously. It is not clear whether the expressions  2*x+x^2  and  sqrt(4*x^3+x^4)  must be in the numerator or in the denominator. I assumed that in the numerator.

restart;
P1:=plots:-shadebetween(0, 1/2*(2*x+x^2)+1/2*sqrt(4*x^3+x^4), x=0..1/13*(5-2*sqrt(3)), color="LightBlue"):
P2:=plots:-shadebetween(1/2*(2*x+x^2)+1/2*sqrt(4*x^3+x^4), 2*x,x=0..1/13*(5-2*sqrt(3)), color="Pink"):
T:=plots:-textplot([[0.1,0.05,c1],[0.1,0.17,c2]], font=[times,bold,16]):
plots:-display(P1,P2,T, view=[0..0.15,0..0.25], labels=[x,z]);

Edit.

restart;
ode:=x*diff(y(x),x)*y(x) = (y(x)^2-9)^(1/2):
ic:=y(exp(4)) = 5:
sol:=dsolve([ode,ic],y(x));
res:=odetest(sol,ode);
solve(res) assuming x>0;

Unfortunately, assuming real  often does not work and a more specific assumption is needed.

Edit.  Below is another solution, in which we first automatically find the domain of the expression  res , and then use this in  assuming :

restart;
ode:=x*diff(y(x),x)*y(x) = (y(x)^2-9)^(1/2):
ic:=y(exp(4)) = 5:
sol:=dsolve([ode,ic],y(x));
res:=odetest(sol,ode);
F:=indets(res, function);
R:=[solve(`and`(seq(`or`(F[i]>0, F[i]<0, F[i]=0), i=1..nops(F))))]; # The domain of res
solve(res) assuming `or`(seq(x in R[i],i=1..nops(R)));

The reason for the error (when calculating the derivative at a specific point) is a premature calculation. Here are 3 correct ways, starting with the best one:

restart;
y := t -> t^2: 
z := D(y);
z(1);

# or

y := t -> t^2:
z:=s->eval(diff(y(t),t), t=s);
z(1);

# or

y := t -> t^2:
z:=unapply(diff(y(t),t), t);
z(1);

restart;

S:=sum(1/x[i],i=1..5);
# or
S:=add(1/x[i],i=1..5);
# or
S:=`+`(seq(1/x[i],i=1..5));

That you ask in the additional question is called a polynomial decomposition (representation of a polynomial as a composition of two polynomials of lesser degree). This is (in general) impossible for arbitrary polynomials  g  and  h . If   f@g=h , then obviously the necessary condition  degree(h)=degree(f)*degree(g)  should be true that is not fulfilled in your example. Maple has the  compoly  command to decompose polynomials. See example below:

restart;
g:=unapply((x+1)^2, x);
f:=unapply(2*x^2+x+3, x);
h:=expand((f@g)(x));
compoly(h, x); 

expand(subs(%[2], %[1]));  # Check

In this example, we see that such a decomposition (if it is possible) is not unique.

Use the  tickmarks  option  for this. See help on  tickmarks  for details .
An example:

plot(x^2, x=-2..2, tickmarks=[0,0]);