We know Frobenius theorem handle pde systems like $\{Xf=0, Yf=0\}$ requiring Lie bracket $[X,Y]\equiv 0 \mod X, Y$ for completely integrability of the system. However, how to handle systems like $\{Xf+Yg=0\}$. For example, how to solve $\frac{\partial u_i}{\partial x_j}=\frac{\partial u_j}{\partial x_i}$ for $i,j=1..n$?

[too long for a comment] Yes, the solution follows precisely from the Poincare' lemma. To see the relationship between Frobenius theorem and the theory of prolongations, let $\pi:\mathbb{R}^2\times\mathbb{R}\to\mathbb{R}^2$ be the trivial bundle, and $x,y$ coordinates on its basis. Let $$ X=a\partial_x+b\partial_y $$ $$ Y=a'\partial_x+b'\partial_y $$ and regard (any) vector field $X$ as a smooth function $\psi_X$ on the first jet extension $J^1(\pi)$, i.e., the space with coordinates $x,y,u,u_x,u_y$: $$ \psi_X=au_x+bu_y. $$ Recall that the function $\psi_X$ is intrinsically defined by the property $$ \psi_X\circ j_1(f)=X(f) $$ for any function $f=f(x,y)$. In this "jet formalism", your system of equations corresponds to the sub manifold $E=\{\psi_X=0,\psi_Y=0\}$ of $J^1(\pi)$. Now we must introduce the first prolongation $E^{(1)}$ of $E$, which is a submanifold of $J^2(\pi)$: to define it, take the total derivatives of the functions $\psi_X$ and $\psi_Y$: $$ D_x(\psi_X)=a_xu_x+au_{xx}+b_xu_y+bu_{xy} $$ $$ D_y(\psi_X)=a_yu_x+au_{xy}+b_yu_y+bu_{yy} $$ $$ D_x(\psi_Y)=a'_xu_x+a'u_{xx}+b'_xu_y+b'u_{xy} $$ $$ D_y(\psi_Y)=a'_yu_x+a'u_{xy}+b'_yu_y+b'u_{yy}. $$ In terms of ideals, $E$ corresponds to the ideal $I=\mathrm{Span}\{\psi_X,\psi_Y\}$ in $C^\infty(J^1(\pi))$, whereas $E^{(1)}$ corresponds to the ideal $I'=\mathrm{Span}\{\psi_X,\psi_Y,D_x(\psi_X),D_y(\psi_X),D_x(\psi_Y),D_y(\psi_Y)\}$ in $C^\infty(J^2(\pi))$. You should also remember that $C^\infty(J^1(\pi))$ is embedded into $C^\infty(J^2(\pi))$: then, the "algebraic formulation" of formal integrability reads $$ I'\cap C^\infty(J^1(\pi)) = I. $$ (Observe that the "$\supseteq$" inclusion always holds). Now, you can easily check that $I'\cap C^\infty(J^1(\pi))=\mathrm{Span}\{\psi_X,\psi_Y,\psi_{[X,Y]}\}$: basically, you show that the only way (besides taking $\psi_X$ or $\psi_Y$ themselves) to get an element of $C^\infty(J^1(\pi))$ out of the generators of $I'$ is by means of the following $C^\infty(J^1(\pi))$linear combination: $$ aD_x(\psi_Y)+bD_y(\psi_Y)a'D_x(\psi_X)b'D_y(\psi_X), $$ which turns out to be exactly $\psi_{[X,Y]}$, i.e., the function associated to the commutator $[X,Y]$. So,"formal integrability" holds if and only if $$ \mathrm{Span}\{\psi_X,\psi_Y,\psi_{[X,Y]}\}=\mathrm{Span}\{\psi_X,\psi_Y\}, $$ which is just a paraphrase of Frobenius' condition. 

