I don't think you have the right way of obtaining the integration factor.
From what I read online, if the factor is 1/P(x,y) then it needs to satisfy
∂(M/P)/∂y=∂(N/P)/∂x

where the inexact equation is M.dx+N.dy=0
Your factor does not seem to have that property.
This is not an area I have ever worked in, so I don't know of a better way; I assumed P was a quadratic in x and y and obtained a solution by sheer algebra.

As a check, if the x2 coefficient is 1 then the constant term in the quadratic is 26.

I don't think you have the right way of obtaining the integration factor.
From what I read online, if the factor is 1/P(x,y) then it needs to satisfy
∂(M/P)/∂y=∂(N/P)/∂x

where the inexact equation is M.dx+N.dy=0
Your factor does not seem to have that property.
This is not an area I have ever worked in, so I don't know of a better way; I assumed P was a quadratic in x and y and obtained a solution by sheer algebra.

As a check, if the x2 coefficient is 1 then the constant term in the quadratic is 26.

I certainly can see a mistake in the OP. You say at some point that the term b is free of x, but if I am not mistaken it is

##b=\frac{-x+y+6}{x^2+y^2+4x+6y}## (another mistake i see is that you seem to say that ##b=\frac{x+y+4}{x^2+y^2+4x+6y}## which is equal to a, I reckon this is a cut and paste mistake ), it clearly is NOT free of x.

I certainly can see a mistake in the OP. You say at some point that the term b is free of x, but if I am not mistaken it is

##b=\frac{-x+y+6}{x^2+y^2+4x+6y}## (another mistake i see is that you seem to say that ##b=\frac{x+y+4}{x^2+y^2+4x+6y}## which is equal to a, I reckon this is a cut and paste mistake ), it clearly is NOT free of x.

Yes I know that, but even if I corrected that mistake,I can't get my results.and I said that because of the method which is to integrate the part of the "b", which is free of x,and no part of "b" is free of x,so here that value is 0.

Interesting...
What I see online is that the integrating factor being a function of x or y only are easy special cases. Indeed, they are the only examples I have found. I note also that the IF in the example satisfies the equation

where 1/P is the IF. But the IF you obtained using the formula in the text does not satisfy this, so I am suspicious that the formula is not general. Certainly it does not yield the IF I got that satisfies the equation above.

Interesting...
What I see online is that the integrating factor being a function of x or y only are easy special cases. Indeed, they are the only examples I have found. I note also that the IF in the example satisfies the equation

where 1/P is the IF. But the IF you obtained using the formula in the text does not satisfy this, so I am suspicious that the formula is not general. Certainly it does not yield the IF I got that satisfies the equation above.

I have just now realized that my equation is not a homogeneous equation so that's why the formula isn't working.I can't believe I missed such a tiny thing. But can you still tell me now how to find the I.F of this non exact,non homogeneous differential equation?

I have just now realized that my equation is not a homogeneous equation so that's why the formula isn't working.I can't believe I missed such a tiny thing. But can you still tell me now how to find the I.F of this non exact,non homogeneous differential equation?

Try my method in post #2. It's a bit of a slog, but it gets there.
Assume the IF is 1/P where P is some quadratic. You can take the x2 coefficient as 1. Find a solution to ∂(M/P)/∂y=∂(N/P)/∂x.
What to do after that I am not sure.

Try my method in post #2. It's a bit of a slog, but it gets there.
Assume the IF is 1/P where P is some quadratic. You can take the x2 coefficient as 1. Find a solution to ∂(M/P)/∂y=∂(N/P)/∂x.
What to do after that I am not sure.

I think you can eliminate the 4 and 6 with the substitutions ##r=x+y+4,~s = -x+y+6## and make a homogeneous equation in ##r## and ##s## out if it. (Homogeneous in the ##f(\lambda x,\lambda y) = \lambda^n f(x,y)## sense).

I'm sure there's a better way, but I embarked on solving the general problem of finding an integration factor for the general first order case: M=Ax+By+C, N=Dy+Ex+F. The IF must be 1/P where P satisfies ∂(M/P)/∂y=∂(N/P)/∂x. I assumed P to be a quadratic in x and y.

I tried differentiating that and could not get back to the original equation.

Edit:
I get ##\frac 12\ln((x-1)^2+(y+5)^2)+\arctan(\frac{x-1}{y+5})##.

The easiest way to make it homogeneous is to adjust each of x and y by a constant: x=X+1, y=Y-5.
(X+Y).dX+(-X+Y).dY=0
Now we can apply the formula in the text to get the IF ##\frac 1{X^2+Y^2}##. This leads easily to the solution ##\frac 12\ln(X^2+Y^2)+\arctan(\frac{X}{Y})=C##.

Well, I don't know if I want to spend much more time on this, but if you choose the constant right so ##f(0,0) = 0## in both your and my solution, they both give the same graph on the 4x4 square. I used implicitplot with Maple:

I submitted the original problem (no change of variable) to Wolfram Alpha. The result appears to be equivalent to what @LCKurtz has here - after I did quite a bit of algebra and some work with verifying equivalent forms of derivatives of arctangent .

I submitted the original problem (no change of variable) to Wolfram Alpha. The result appears to be equivalent to what @LCKurtz has here - after I did quite a bit of algebra and some work with verifying equivalent forms of derivatives of arctangent .