We consider (3+1)-dimensional second-order evolutionary PDE5 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 the symplectic Monge-Ampere form and determine their Lagrangians. We develop a calculus for transforming the symmetry condition to the "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 x 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. (C) 2018 Elsevier B.V. All rights reserved.