a:=fsolve(1/FINT2);
#                       a := 0.02355400734
e:=1e-6:
int(FINT2, x = 0 .. a-e, numeric);
#                         -0.7603068518
int(FINT2, x = a+e .. 1, numeric);
#                          1.155767949
%+%%;
#                          0.3954610972

You may use the operator D; just fill its remember table

D(y):=1: D(x):=x*y: D(z):=z^2:
D(x*y + y*z);
D(1/2*x*y^2 + 3*y^3*z + a); # you may want D(a):=0;

y(x) = 0 is obtained for _C1 --> infinity.
This is documented (?dsolve,details).

What you write here is simply wrong.

ode:=diff(y(x), x)/y(x);
sol:=y(x)=0;

algsubs(sol,ode);   #this gives ZERO. It should be 1

algsubs should return an error (division by 0). 1 is out of the question.

ode:=diff(y(x),x)^2+2*x*diff(y(x),x)/y(x)-1 = 0;
sol:=y(x)=0;
odetest(sol,ode);
         0

is also wrong (division by 0 too). y(x)=0 is not a solution (odetest checks the solution for a transformed ode, which is not equivalent with the original).

You say: Which is 0/0 but this is 1 in the limit. This is nonsense.

You have changed the problem, now you want  f(t*x,t^p*y) = t^(p-1)*f(x,y).
In this case the procedure is even simpler:

Q := (F,x,y) -> simplify( (F + diff(F,x)*x)/(F - diff(F,y)*y) ):

Let's compare with the previous P:

P := (F,x,y)-> simplify(((-x*diff(F,x,x)-diff(F,x))*F+diff(F,x)^2*x)/(y*(F*diff(F,x,y)-(diff(F,x))*(diff(F,y))))):

Q(-2*x/3+sqrt(x^2+3*y)/3, x,y ),   P(-2*x/3+sqrt(x^2+3*y)/3, x,y );
#                              2, 2

Q(3*sqrt(x*y), x,y),   P(3*sqrt(x*y), x,y);
#                             3, -1

P := (F,x,y)-> simplify(((-x*diff(F,x,x)-diff(F,x))*F+diff(F,x)^2*x)/(y*(F*diff(F,x,y)-(diff(F,x))*(diff(F,y))))):
P(-2*x*(1/3)+(1/3)*sqrt(x^2+3*y), x,y);
#                              2
P(x*y^3-x^2/y^3, x,y);
#                               1/6

I suspect that n was previously assigned to 0.
Execute a restart or n := 'n';

f := sqrt(x):    # assumed to be unknown
u := sprintf("%Zm", f):

sscanf(u, "%m");

==> [sqrt(x)]

(Actually I don't know why "%Zm" was used instead of "%m")

Your expression v depends on x and y. It has not an elementary antiderivative int(v. x) but you can compute easily a definite double integral or a simple one (giving a numerical value to y):

restart;
# https://en.wikipedia.org/wiki/Hilbert_curve
rot:=proc (n, X, Y, rx, ry)
  local x:=X, y:=Y;
  if ry = 0 then
    if rx = 1 then x := n-1 - x; y := n-1 - y fi;
  return y,x;
  fi;
  x,y
end:
xy2d:=proc(n, X, Y)
  local rx, ry, s:=iquo(n,2), d:=0, x:=X,y:=Y;
  while s>0 do
    rx := `if`(Bits:-And(x,s,bits=n) > 0,1,0);
    ry := `if`(Bits:-And(y,s,bits=n) > 0,1,0);
    d += s * s * Bits:-Xor(3 * rx,ry, bits=n) ;
    x,y := rot(n, x, y, rx, ry);
    s:=iquo(s,2)
  od;
  d;
end:
d2xy:=proc(n, d)
  local rx, ry, s:=1, t:=d, x:=0, y:=0;
  while s<n   do  
    rx := irem(iquo(t,2),2);
    ry:=irem(irem(t,2)+irem(rx,2),2);
    x,y:=rot(s, x, y, rx, ry);
    x += s * rx;
    y += s * ry;
    t :=iquo(t,4);
    s *=2;
  od;
  x,y
end:

See at  the application center https://www.maplesoft.com/applications/view.aspx?SID=154350

They are almost equivalent; actually sol2 is better because it is valid for x=0 too.
To answer why solve prefers sol1 would imply to see its code, but probably solve used the substitution x*y = Y (which is natural).

It seems that you actually have normalized barycentric coordinates.
Anyway,

delta:=proc(ex, X,Y,Z)
local c:=map2(coeff, ex/(X*Y*Z),[1/X,1/Y,1/Z]);
(c[2]+c[3]-c[1])^2-4*c[2]*c[3];
end:

delta(2*Y*Z+3*Z*X+4*X*Y, X,Y,Z);
#                              -23

I'd suggest something like this:

restart;
res:=-2*x*csgn(1/(x+1))/(x+1)^3+2*x/(x+1)^3 + csgn((x^2-4))-1;
u:=op~(indets(res,csgn(anything)));
A:=solve~(u) union rhs~(op~(singular~(u)))  union {-infinity, infinity};
A:=sort([select(type,A, realcons)[]]);
for i to nops(A)-1 do
   A[i] .. A[i+1], simplify(res) assuming x>A[i],x<A[i+1]
od;

                                      4 x   
                   -infinity .. -2, --------
                                           3
                                    (x + 1) 
                              3      2          
                          -2 x  - 6 x  - 2 x - 2
                -2 .. -1, ----------------------
                                        3       
                                 (x + 1)        
                          -1 .. 2, -2
                        2 .. infinity, 0

There are (in general) many possibilities to write a polynomial as a sum of squares.
For your example, here are some such results, including the desired one.

restart;
ex:=expand((x - y)^2 + (y - z)^2 + (z - x)^2):
ex1:=add((a[i]*x+b[i]*y+c[i]*z)^2, i=1..3):
zzz:=coeffs(expand(ex-ex1),[x,y,z]):
sols:=solve([zzz, a[1]=1,a[2]=1,b[3]=1/2], explicit):
'ex'=eval(ex1, sols[1]);
sols:=solve([zzz, a[1]=1,a[2]=1/2,b[3]=1/2], explicit):
'ex'=eval(ex1, sols[1]);
sols:=solve([zzz, a[1]=1,b[3]=0,a[2]=0], explicit):
'ex'=eval(ex1, sols[1]);