@vv It seems that the answer you have given identifies the homogeneous part of the polynomial solution post processedly. But I want the program to look for the homogeneous(with eliminated terms) solutions beforehand so that the computing time is reduced when I want to calculate the polynomial solutions of higher degree.

@Carl Love Thanks for the reply.

@Carl Love For example

PDE := 2*((x1 - x2)^2 + (y1 - y2)^2)*(u1*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), x1) + v1*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), y1) + u2*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), x2) + v2*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), y2)) - ((u1 - v1)^2 + (u2 - v2)^2)*((x1 - x2)*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), u1) + (y1 - y2)*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), v1) - (x1 - x2)*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), u2) - (y1 - y2)*diff(f(x1, y1, x2, y2, u1, v1, u2, v2), v2)) = 0

PolynomialSolutions(PDE, degree = 1)  gives

{f(x1, y1, x2, y2, u1, v1, u2, v2) = (v1 + v2)*_C2 + (u1 + u2)*_C3 + _C1}

PolynomialSolutions(PDE, degree = 2) gives

{f(x1, y1, x2, y2, u1, v1, u2, v2) =_C1 + (v1 + v2)*_C2 + (u1 + u2)*_C3 + (v1 + v2)^2*_C4 + (v1 + v2)*(u1 + u2)*_C5 + (u1 + u2)^2*_C6 + (u1*y2 + u2*y1 - v1*x2 - v2*x1)*_C7 + (u1*y1 + u2*y2 - v1*x1 - v2*x2)*_C8}

But from degree 1 solution we already know that (u1 + u2) , (v1 + v2)  and their functions are solutions of the pde. So we dont need the terms (u1 + u2) , (v1 + v2) , (u1 + u2)^2 , (v1 + v2)^2 , (v1 + v2)*(u1 + u2) in degree 2 solution because we already know that they are solutions. We only need the terms (u1*y2 + u2*y1 - v1*x2 - v2*x1) and (u1*y1 + u2*y2 - v1*x1 - v2*x2). So, how do we find only these new solutions without having to calculate the terms which we already know are solutions?