If we can ignore singular behavior ([tex]\vec{N}\neq 0[/tex], [tex]b(x)-u(x) \neq 0[/tex], etc), one way to try to solve
[tex]\nabla u = (b-u)\vec{N}[/tex]
is the following. Define a new orthonormal basis
[tex]\vec{v} = \frac{1}{\sum_i N_i^2} (N_1, \ldots, N_n),[/tex]
[tex]\vec{v}^{(1)}_\perp = \frac{1}{N_1^2+N_2^2} (N_2, -N_1,0,\ldots,0).[/tex]
The other [tex]\vec{v}^{(i)}_\perp[/tex] can be constructed by Gram-Schmidt. So these are all an orthonormal set.
We then note that
[tex]\vec{v} \cdot \nabla u = (b-u)[/tex]
[tex]\vec{v}^{(i)}_\perp \cdot \nabla u = 0.[/tex]
If we can define new coordinates [tex](X,Y_i)[/tex] by inverting
[tex]\frac{\partial}{\partial X} = \vec{v} \cdot \nabla_x,~~\frac{\partial}{\partial Y_i} = \vec{v}^{(i)}_\perp \cdot \nabla_x,[/tex]
then we obtain a first-order ordinary differential equation for [tex]u(X)[/tex]:
[tex]\frac{\partial u}{\partial X} = b- u, ~~ \frac{\partial u}{\partial Y^i} =0.[/tex]
I suspect that one of the requirements for these new coordinates to exist is that [tex]N = \nabla\phi[/tex] for some scalar [tex]\phi[/tex] so that 2nd derivatives match.