In this paper, we study a new numerical method of order 4 in space and time based on half-step discretization for the solution of three-space dimensional quasi-linear hyperbolic equation of the form utt=A(x,y,z,t,u)uxx+B(x,y,z,t,u)uyy+C(x,y,z,t,u)uzz+f(x,y,z,t,u,ux,uy,uz,ut), A>0, B>0, C>0, 0<x, y, z<1, t>0 subject to prescribed appropriate initial and Dirichlet boundary conditions. The scheme is compact and requires nine evaluations of the function f. For the derivation of the method, we use two inner half-step points each in x-, y-, z- and t-directions. The proposed method is directly applicable to three-dimensional (3D) hyperbolic equations with singular coefficients, which is main attraction of our work. We do not require extra grid points for computation. The proposed method when applied to 3D damped wave equation is shown to be unconditionally stable. Operator splitting method is used to solve 3D linear hyperbolic equations. Few benchmark problems are solved and numerical results are provided to support the theory which is discussed in this paper.