Your terminology is a little off, the set $I_c$ is not called a "level function" it is called a "level set" or, as you say in your title, a "level surface". The intuition is that it consists of all points at the same "level", where different values of $c$ define different levels.
In general there's not much you can say about the topology of $I_c$ to justify calling it a level surface. However, if $c$ is a regular value of the function then you can apply the submersion theorem, which is basically a version of the implicit function theorem, to prove that $I_c$ is a submanifold of dimension $n-1$. In this case one might say that $I_c$ is a "level manifold" of the function, in the case $n=3$ it would truly be a "level surface", but even in higher dimensions the terminology is abused and $I_c$ is called a "level surface", or you could simply call it a level hypersurface without any abuse of terminology.
If all values of $c$ are regular values, and if $\lambda$ is surjective, then yes, the set of level manifolds forms a foliation, and it is a foliation of dimension $n-1$, where the dimension refers to the dimensions of the individual leaves of the foliation, i.e. the individual level manifolds.
But if some values of $c$ are not regular values, or if $\lambda$ is not surjective, then generally speaking you cannot call the collection of level sets a "foliation" (except in highly special circumstances).
By the way, there are plenty of good books which cover these topics in detail. I like the first volume of Spivak's "Differential Geometry" for these topics in your post.
Let $x\in \mathbb{R}^n$ and $U$ be an open neighbourhood of $x.$ Let $\Delta$ be a smooth $k$-dimensional distribution on $\mathbb{R}^n.$ When $\Delta$ is involutive, we have, by Frobenius' theorem, that it is locally completely integrable. That means that there exists a coordinate transformation where the immersed submanifolds tangent to $\Delta$ are "flattened" in the new coordinates.
Let us use that change of coordinates. Let the coordinate transformation be denoted $\Phi: U\to V.$ Define our new set of coordinates $$\begin{pmatrix}z_1(x) \\ \vdots \\z_n(x)\end{pmatrix} = z(x) = \Phi(x).$$ The sets tangent to $\Delta$ are immersed submanifolds given in the new coordinates $z$ by fixing the $n-k$ functions $z_{k+1}(x), \dots, z_n(x)$ to any constant. These are your $\lambda$ functions.
It helps to move to these new coordinates $z$ where $\Delta$ is flattened. Notice that the $\Phi$-related distribution of $\Delta$ is generated by the vector fields
$$\partial_{z_1},\dots,\partial_{z_k}.$$
Let us call this distribution (defined on the open set $V$) $\bar{\Delta}.$ Let us also denote the $\Phi$-related vector field of $f$ as $\bar{f}.$
All of this discussion has ignored the important property connecting $f$ and $\Delta.$ Now let us talk about that.
Since $\Delta$ is involutive we also have $\bar{\Delta}$ is involutive. Further, if $[f, \Delta] \subseteq \Delta$ we have that $[\bar{f}, \bar{\Delta}] \subseteq \bar{\Delta}$
Recognize that since $\bar{\Delta}$ is generated by the constant, standard vector fields $\partial_{z_1}, \dots, \partial_{z_k}$ we can say
$$
\begin{aligned}\\
[\bar{f}, \partial_{z_1}] &= \sum_{\ell=1}^{k} c_{1,\ell} \partial_{z_\ell}\\
&~\vdots\\
[\bar{f}, \partial_{z_k}] &= \sum_{\ell=1}^{k} c_{k,\ell} \partial_{z_\ell}
\end{aligned}$$
where $c_{i,j}$ are smooth functions on $V.$ At this point if you write $\bar{f}$ as a smooth function combination of the constant vector fields $\partial_{z_1},\ldots,\partial_{z_n}$ and combine with the above equation, what can you say about the coefficients that multiply the vector fields $\partial_{z_{k+1}},\dots, \partial_{z_n}$? Direct computation should verify that those coefficients cannot be functions of $z_{1}$ through til $z_k.$
Best Answer
One way to solve such systems is given in the proof of the Frobenius theorem. The proof of sufficiency is constructive; it provides a way to construct a solution. There is, however, an easier way.
A system consisting of one equation can be solved using the method of characteristics. Consider a system consisting of one (first) equation: $$\tag{1} \frac{\partial \lambda }{\partial x_1}x_1-\frac{\partial \lambda }{\partial x_2}x_2=0 $$ Note that on the left side of (1) there is nothing other than the derivative of the function $\lambda(x_1,x_2,x_3)$ along the trajectories of an ode system $$\tag{2} \frac{dx_1}{dt}= x_1\qquad \frac{dx_2}{dt}= -x_2\qquad \frac{dx_3}{dt}= 0 $$ and $\lambda(x_1,x_2,x_3)$ satisfies (1) iff it is a first integral of (2). Let's find the first integrals of(2). $$ dt=\frac{dx_1}{x_1}=\frac{dx_2}{-x_2},\qquad dx_3=0 $$ gives us two first integrals $$ \phi_1(x_1,x_2,x_3)= x_1 x_2,\qquad \phi_2(x_1,x_2,x_3)= x_3. $$ The general solution of (1) consists of all first integrals of (2), i.e. $$ \lambda(x_1,x_2,x_3)=F(\phi_1(x_1,x_2,x_3),\phi_2(x_1,x_2,x_3))= F(x_1x_2,x_3) $$ where $F(\cdot,\cdot)$ is any function of the required order of continuity.
Substitute the obtained general solution into the second equation $$\tag{3} \frac{\partial \lambda }{\partial x_1}x_1-\frac{\partial \lambda }{\partial x_3}x_3=0. $$ Let $\lambda(x_1,x_2,x_3)= F(u,v)$, $u=x_1x_2$, $v=x_3$. Using the chain rule one obtains $$ \frac{\partial\lambda}{\partial x_1}= \frac{\partial F}{\partial u}\frac{\partial u}{\partial x_1}+ \frac{\partial F}{\partial v}\frac{\partial v}{\partial x_1}=\frac{\partial F}{\partial u} x_2 $$ ($\frac{\partial\lambda}{\partial x_2}$ is not needed) $$ \frac{\partial\lambda}{\partial x_3}= \frac{\partial F}{\partial u}\frac{\partial u}{\partial x_3}+ \frac{\partial F}{\partial v}\frac{\partial v}{\partial x_3}=\frac{\partial F}{\partial v}. $$ Substitution into (3) gives us $$ \frac{\partial F}{\partial u} x_2x_1-\frac{\partial F}{\partial v}x_3=0 $$ $$\tag{4} \frac{\partial F}{\partial u} u-\frac{\partial F}{\partial v}v=0 $$ The equation (4) can also be solved using the method of characteristics: $$ \frac{du}{u}=\frac{dv}{-v} $$ implies that $$ F(u,v)= G(uv), $$ where $G(\cdot)$ is any function of the required order of continuity. So, the general solution to the original system is $\lambda(x_1,x_2,x_3)= G(x_1x_2x_3)$.