Using the equivalences in the comments, we can also write the following partial derivatives with some manipulation:
r=∂p/∂x, s=∂p/∂y and s=∂q/∂x, t=∂q/∂y from the definitions of p=∂z/∂x and q=∂z/∂y.
p and q are functions of x and y only, so only partial derivatives wrt x and y are required.
And we can relate partial to full derivatives (for p(x,y) and q(x,y)):
dz=(∂z/∂x)dx+(∂z/∂y)dy=pdx+qdy; similarly:
dp=rdx+sdy, dq=sdx+tdy. From these we get: r=(dp-sdy)/dx and t=(dq-sdx)/dy.
The given PDE needs to be written in the form: Rr+Ss+Tt=V:
-xyr+(x2-y2)s+xyt=2x2-2y2+py-qx, so:
R=-xy, S=x2-y2, T=xy, V=2x2-2y2+py-qx.
Rr+Ss+Tt=R(dp-sdy)/dx+Ss+T(dq-sdx)/dy=V,
Rdy(dp-sdy)+Ssdxdy+Tdx(dq-sdx)=Vdxdy,
Rdydp-sRdy2+sSdxdy+Tdxdq-sTdx2-Vdxdy=0,
Rdydp+Tdxdq-Vdxdy=s(Rdy2-Sdxdy+Tdx2).
One possible solution to these equations is:
Rdydp+Tdxdq-Vdxdy=0 and Rdy2-Sdxdy+Tdx2=0 (Monge's subsidiary equations).
That is:
-xydydp+xydxdq-(2x2-2y2+py-qx)dxdy=0 and:
-xy(dy)2-(x2-y2)dxdy+xy(dx)2=0,
(dy)2+(x/y-y/x)dxdy-(dx)2=0.
Consider (dy-m1dx)(dy-m2dx)=(dy)2-(m1+m2)dxdy+m1m2(dx)2=0.
So m1+m2=y/x-x/y and m1m2=-1; m1=y/x, m2=-x/y, making m1m2=-1.
dy-(ydx/x)=0, xdy=ydx; or, dy+xdx/y=0, ydy=-xdx. We have two auxiliary integrals:
(1) dy/y-dx/x=0; integrating: ln|y|-ln|x|=k, ln|y/x|=k, or y/x=ek, which we can write y/x=a, a constant, so y=ax and dy=adx.
(2) ydy+xdx=0; integrating: ½y2+½x2=c, a positive constant (because the sum of the squares of two real numbers is always positive), which can be written x2+y2=b2, where c=b2.
Therefore, ydy=-xdx, and y=√(b2-x2), so dy=-xdx/y=-xdx/√(b2-x2).
We can use these to make substitutions in the other subsidiary equation.
(1)⇒(1i) -x(ax)(adx)dp+x(ax)dxdq-(2x2-2(ax)2+p(ax)-qx)dx(adx)=0, (substituting y=ax and dy=adx)
-a2x2dxdp+ax2dxdq-a(2x2-2a2x2+apx-qx)(dx)2=0,
-a2x2dxdp+ax2dxdq-2ax2(dx)2+2a3x2(dx)2-a2px(dx)2+aqx(dx)2=0,
axdx(-axdp+xdq-2xdx+2a2xdx-apdx+qdx)=0
-axdp+xdq-2xdx+2a2xdx-apdx+qdx=0,
-a(xdp+pdx)+(xdq+qdx)+2xdx(a2-1)=0.
This is integrable: -apx+qx+(a2-1)x2=c, a constant.
We have a=y/x and c=-apx+qx+(a2-1)x2, so, if c=f1(a), where f1 is an arbitrary function:
-apx+qx+(a2-1)x2=f1(y/x).
(2)⇒(2i) -xydydp+xydxdq-(2x2-2y2+py-qx)dxdy=0 is the subsidiary equation.
Substitute y=√(b2-x2) and dy=-xdx/√(b2-x2) (also ydy=-xdx).
-x(-xdx)dp+x√(b2-x2)dxdq-(2x2-2(b2-x2)+p√(b2-x2)-qx)dx(-xdx/√(b2-x2))=0.
We can divide through by xdx:
xdp+√(b2-x2)dq+(4x2-2b2+p√(b2-x2)-qx)dx/√(b2-x2)=0,
xdp+√(b2-x2)dq+(4x2-2b2)dx/√(b2-x2)+pdx-qxdx/√(b2-x2)=0,
d(px)+d(q√(b2-x2))+(4x2-2b2)dx/√(b2-x2)=0.
Integrating: px+q√(b2-x2)+∫[(4x2-2b2)/√(b2-x2)]dx=c, a constant.
Let F(x)=(4x2-2b2)dx/√(b2-x2), then c=px+q√(b2-x2)+F(x). (F(x) is integrable, but to save space the integration is not shown here.)
We have b2=x2+y2 and c=px+q√(b2-x2)+F(x). If c=f2(b2) then:
px+q√(b2-x2)+F(x)=f2(x2+y2). But -apx+qx+(a2-1)x2=f1(y/x), q=(apx+f1(y/x)-(a2-1)x2)/x from (1i).
px+(apx+f1(y/x)-(a2-1)x2)/x)√(b2-x2)+F(x)=f2(x2+y2),
(px2+apx+f1(y/x)-(a2-1)x2)√(b2-x2)+xF(x)=xf2(x2+y2),
p(x2+ax)√(b2-x2)=xf2(x2+y2)-xF(x)-(f1(y/x)-(a2-1)x2)√(b2-x2),
p=(xf2(x2+y2)-xF(x)-(f1(y/x)-(a2-1)x2)√(b2-x2))/((x2+ax)√(b2-x2)).
dz=pdx+qdy, so p and q can be substituted (!) to get dz in terms of x and y.
The above calculations need to be verified and corrected if necessary...