We present a new two-level implicit difference method of $O(k^2 + kh^2 + h^4)$ for approximating the three space dimensional non-linear parabolic differential equation $u_{xx} +u_{yy} +u_{zz} = f(x, y, z, t, u, u_x, u_y, u_z, u_t)$, $0<x, y, z<1$, $t>0$ subject to appropriate initial and Dirichlet boundary conditions, where $h>0$ and $k>0$ are mesh sizes in space and time directions, respectively. In addition, we also propose some new two-level explicit stable methods of $O(kh^2+h^4)$ for the estimates of ($∂u/∂n$). When grid lines are parallel to $x-$, $y-$ and $z-$ coordinate axes, then ($∂u/∂n$) at an internal grid point becomes ($∂u/∂x$), ($∂u/∂y$) and ($∂u/∂z$), respectively. In all cases, we require only 19-spatial grid points and a single computational cell. The proposed methods are directly applicable to singular problems and we do not require any special technique to handle singular problems. We also discuss operator splitting method for solving linear parabolic equation. This method permits multiple use of the one-dimensional tri-diagonal solver. It is shown that the operator splitting method is unconditionally stable. Numerical tests are conducted which demonstrate the accuracy and effectiveness of the methods developed.