  We consider (3+1)-dimensional second-order evolutionary PDEs where the
unknown $u$ enters only in the form of the 2nd-order partial derivatives. For
such equations which possess a Lagrangian, we show that all of them have a
symplectic Monge--Amp\`ere form and determine their Lagrangians. We develop a
calculus for transforming the symmetry condition to a "skew-factorized" form
from which we immediately extract Lax pairs and recursion relations for
symmetries, thus showing that all such equations are integrable in the
traditional sense. We convert these equations together with their Lagrangians
to a two-component form and obtain recursion operators in a $2\times 2$ matrix
form. We transform our equations from Lagrangian to Hamiltonian form by using
the Dirac's theory of constraints. Composing recursion operators with the
Hamiltonian operators we obtain the second Hamiltonian form of our systems,
thus showing that they are bi-Hamiltonian systems integrable in the sense of
Magri. By this approach, we obtain five new bi-Hamiltonian multi-parameter
systems in (3+1) dimensions.
