The numerical and analytic solutions of the mixed problem for multidimensional fractional hyperbolic partial differential equations with the Neumann condition are presented. The stable difference scheme for the numerical solution of the mixed problem for the multidimensional fractional hyperbolic equation with the Neumann condition is presented. Stability estimates for the solution of this difference scheme and for the first- and second-order difference derivatives are obtained. A procedure of modified Gauss elimination method is used for solving this difference scheme in the case of one-dimensional fractional hyperbolic partial differential equations. He's variational iteration method is applied. The comparison of these methods is presented.