Solving by Lagrange's auxiliary equations $ \frac {dx}{z} =\frac {dy}{z} =\frac {(3z+1)(dz)}{x+y} $ By taking 1st and 2nd equations 1st solution is $$ x=y+c_{0}$$
but i am getting two different solutions slightly differing :
Solution 1 :  using lagrange multiplier 1,1,0 we get
$$ \frac{dx+dy}{2z}=\frac{(3z+1)dz}{x+y}$$
$$ \frac {(x+y)d(x+y)}{2z}= {(3z+1)(dz)}$$
$$ \frac {d(x+y)^{2}}{2} = (2z)(3z+1)(dz) =(6z^{2}+2z)dz$$
$$ \frac {(x+y)^2}{2} = {2z^3}+(z^2)+c_{1} $$
$$ (x+y)^{2}= 4z^{3}+2z^{2}+c_{1'}    \to(1)$$
Solution 2 : taking x,y,-z as multipliers
$$ \frac {xdx+ydy-(zdz)(3z+1)}{xz+yz-z(x+y)} =\frac {xdx+ydy-(zdz)(3z+1)}{0}$$
$$ xdx+ydy= (3z+1)(zdz)=(3z^{2}+z)dz$$
$$ \frac{x^2}{2} +\frac{y^2}{2}= z^{3}+ \frac{z^2}{2}+c_{2} $$
$$ x^{2}+y^{2}=2z^{3}+z^{2}+c_{2'}$$
$$ 2(x^{2}+y^{2})=4z^{3}+2z^{2}+c_{2''}                \to(2)$$
here the first solution $(1)$ is different  and the term of $$xy$$ is extra .. i believe i have done mistake in my first solution but i am unable to find out where?
solution $(1)$ is my solution and solution $(2)$ is the solution from the the institute

There is no mistake in your calculus.
Both :
$$ (x+y)^{2}= 4z^{3}+2z^{2}+c_{1'}    \to(1)$$
$$ 2(x^{2}+y^{2})=4z^{3}+2z^{2}+c_{2''}                \to(2)$$
are correct, each one on different characteristic curves (because different constants of integrations).
This is not contradictory (on the characteristic curves, but not elsewhere) because
$$4z^3+2z^2=(x+y)^2-c_{1''}=2(x^2+y^2)-c_{2''}$$
$$2(x^2+y^2)-(x+y)^2=c_{2''}-c_{1''}$$
$$x^2+y^2-2xy=c_{2''}-c_{1''}$$
$$(x-y)^2=c_{2''}-c_{1''}$$
Since $x-y=c_0$
$$c_0^2=c_{2''}-c_{1''}$$
This denotes the relationship between the constants of integrations as expected.
