Suppose in $\mathbf{R}^N$ we have a one-form field,
$ \theta = \sum_{i=1}^N \theta_i d x_i $.
The Pfaff integrability condition is
$d \theta \wedge \theta = 0$.
Is it possible to give an elementary, elegant proof of this theorem without using differential geometry? I mean, how to prove it for first- or second-year students?
You can do it using the theory of (ordinary) differential equations.
Suppose you have the Pfaff form $\theta = \sum \theta_i d x_i$. It gives at every point the normal to the tangent of a (hypothetical) hypersurface $F(x_1, \ldots, x_{n+1}) = \text{const}$, and we want to show that ( given the integrability conditions) such a hypersurface exists through every point ( locally).
Now, try to find the hypersurface as a graph of a function $$x_{n+1} = \phi(x_1, \ldots, x_{n})$$
The condition $\sum \theta_i d x_i= 0$ can also be written as $$d x_{n+1} = \sum_{i=1}^{n} - \frac{\theta_i}{\theta_{n+1}} d x_i$$
which is $d$ of the equality $x_{n+1} = \phi(x_1, \ldots , x_{n})$. Now, to find $\phi$ we have to solve the (Frobenius ) system for the function $\phi(x_1, \ldots, x_n)$
$$\frac{\partial \phi}{\partial x_i} = -\frac{\theta_i(x_1, \ldots, x_{n-1}, \phi(x_1, \ldots, x_{n}))}{\theta_{n+1}(x_1, \ldots, x_{n}, \phi(x_1, \ldots, x_{n}))}= a_i(x_1, \ldots, x_{n}, \phi(x_1, \ldots, x_{n}))$$ for $i=1, \ldots, n$
We've seen such systems before, where the RHS $a_i$, does not depend on $\phi$. The integrability condition is $\frac{\partial a_i}{\partial x_j} = \frac{\partial a_j}{\partial x_i}$, coming from $\frac{\partial^2 \phi}{\partial x_i \partial x_j} = \frac{\partial^2 \phi}{\partial x_j \partial x_i}$. In our case the same condition translates into some conditions for $a_i$, or $\theta_1$, $\ldots$, $\theta_{n+1}$ which are equivalent to $d \theta \wedge \theta = 0$. We skip the details for now.
Note: we can have systems of Pfaff equations and again the problem is if we have a solution of maximal possible dimension ( a subvariety whose dimension is maximal ( if the Pfaff system consists of $d$ forms, we want the variety of codimension $d$). The Pfaff condition expressed with differential forms is more elegant than the Frobenius one, and is also seen to be invariant ( I wonder if it originates with E. Cartan).
Let's sketch an example from the book of Forsyth, Theory of Differential Equations, vol I.
$$\theta = (y+z)(z+u)(u+y) d x + (x+z)(z+u)(u+x) d y + (x+y)(y+u)(u+x) d z + (x+y)(y+z)(z+x) du $$ (symmetric in $x$, $y$, $z$, $u$, so we expect as solution a symmetric function $F(x,y,z,y) = \text{const}$). But now we break some symmetry and consider the Frobenius system in $u = u(x,y,z)$
$$\frac{\partial u}{\partial x} = - \frac{(u+z)(u+y)}{(x+y)(x+z)} \\ \frac{\partial u}{\partial y} = - \frac{(u+x)(u+z)}{(y+x)(y+z)}\\ \frac{\partial u}{\partial z} = - \frac{(u+x)(u+y)}{(z+x)(z+y)}$$
With Mathematica we get a solution
$$u= - \frac{x y z - c(x+y+z)}{x y + y z + x z - c}$$
or, equivalently,
$$F(x,y,z,u) \colon =\frac{u x y + u x z + u y z + x y z}{x+y + z} = c = \text{const}$$
So what we see is $\theta$ is proportional to $d F(x,y,z,u)$
Note: if the condition $d \theta \wedge \theta = 0$ were not satisfied, the Frobenius system would have no solution. So one could go ahead and try to solve it anyways.