Usually, planar and spatial problems of the theory of elasticity in stresses are solved using the Airy or Maxwell stress function, respectively. Solving a plane problem usually comes down to solving a biharmonic equation. In the spatial case, it can be solved by the Castigliano variational method. This work is devoted to the mathematical and numerical solution of the spatial problem of elasticity theory directly concerning stresses. In this case, the boundary value problem consists of three equilibrium equations and three Beltrami-Michell equations with the corresponding boundary conditions. At the same time, equilibrium equations are considered as missing boundary conditions. Symmetric finite-difference equations are constructed, and the problems of equilibrium of a parallelepiped under the action of domed and uniformly distributed loads are solved numerically. The reliability of the results obtained is substantiated by comparing the numerical results with the wellknown Filonenko-Borodich solution and with the results of solving the boundary value problem of the equilibrium of the parallelepiped with respect to displacements.