Is every constraint involving only two coordinates integrable? There is a footnote on Goldstein's Classical Mechanics (3rd ed., page 15) which says the following:

In principle, an integrating factor can always be found for a first-order differential equation of constraint in systems involving only two coordinates and such constraints are therefore holonomic.

I am trying to show this but I can only do that for linear constraints. In that case the constraint can be written as
$$\frac{dy}{dx}+a(x)y=b(x),\tag1$$
and after multiplying both sides by an integrating factor $e^{\left(\int a(x)dx\right)}$, one readily gets the integrable equation
$$d\left(e^{\left(\int a(x)dx\right)}y\right)=b(x)e^{\left(\int a(x)dx\right)}dx.$$
The point is I cannot write constraints such as 
$$dy+\left[a(x)y^2+b(x)\right]dx=0,\tag2$$
in the form of $(1)$. Is there something missing in Goldstein's claim? If not, how to prove it for non-linear constraints?
 A: (Assuming the usual continuity and nonsingularity conditions) A 1st order differential constraint of two variables $x$ and $y$ can be written as $$\begin{matrix}\omega = M(x,y)dx + N(x,y)dy=0& [1]
\end{matrix}$$ 
Rearrange [1] to $$\begin{matrix}\frac{dy}{dx} = -\frac{M(x,y)}{N(x,y)} & [2]
\end{matrix}$$ that leads to the form of a standard ordinary 1st order differential equation in one dependent variable $y$ and one independent variable $x$, and one that has a solution with a given initial condition $x_0,y_0$ in the implicit form of $F(x,y)=k$ for some $k$ that depends on the initial condition. This implicit equation is a solution to [1] and then is equivalent to $dF=0$ or in coordinates $$\begin{matrix}\frac{\partial F}{\partial x}dx+\frac{\partial F}{\partial y}dy=0 & [3]\end{matrix}$$
But [1] and [3] can exist simultaneously iff  $\frac{\frac{\partial F}{\partial x}}{M(x,y)}=\frac{\frac{\partial F}{\partial y}}{N(x,y)}=\lambda(x,y)$ for some $\lambda$ that is just the sought for integrating multiplier, i.e. $$\begin {matrix}dF=\lambda\omega & [4]\end {matrix}$$
A: Comments to the post (v2): 


*

*Be aware that Goldstein later in eq. (2.20') discusses more systematically semi-holonomic constraints that are allowed to depend on time $t$. However on page 15 Goldstein assumes implicitly that there is no explicit $t$-dependence, i.e. there are only 2 generalized coordinates, say $x$ and $y$, without time $t$. The semi-holonomic constraint next kills 1 d.o.f.

*Goldstein is a physics (rather than a math) textbook. The reader is expected to figure out the mathematimatical shortcomings by himself! E.g. there are implicit regularity condition (such as, e.g., differentiability). See also this and this related Phys.SE posts.

*In particular, there can be global obstructions. What Goldstein is trying to convey is the fact that an inexact differential 
$$\omega~=~f(x,y)\mathrm{d}x+g(x,y)\mathrm{d}y\tag{A}$$ 
on a 2-dimensional manifold $M$ (for a point $p\in M$ with $\omega_p\neq 0$ non-vanishing) has an integrating factor $\lambda$ (in a sufficiently small open neighborhood $U\ni p$) that makes the one-form $\left. \lambda\omega\right|_U$ exact. The condition 
$$\left(f\frac{\partial }{\partial y}- g\frac{\partial }{\partial x}\right)\ln\lambda~=~ \frac{\partial g}{\partial x}-\frac{\partial f}{\partial y}\tag{B}$$
for $\ln\lambda$ is a linear first-order PDE in 2 variables, which one may show has local solutions. An exact differential corresponds to a holonomic constraint. 
