Alternative discretization and solution procedures are developed for the 1-D dual phase-lag (DPL) equation, a partial differential equation for very short time, microscale heat transfer obtained from a delay partial differential equation that is transformed to the usual non-delay form via Taylor expansions with respect to each of the two time delays. Then in contrast to the usual practice of decomposing this equation into a system of two equations, we utilize this formulation directly. Truncation error analysis is performed to show consistency and first-order temporal accuracy of the discretized 1-D DPL equation, and it is shown by von Neumann stability analysis and numerical results that the proposed numerical technique is unconditionally stable. The overall approach is then extended to three dimensions via Douglas-Gunn time-splitting, and a simple argument for stability is given for the time-split formulation. Based on a straightforward arithmetic operation count, qualitative comparisons are made with explicit and iterative methods with the expected result that the current approach is generally significantly more efficient, and this is demonstrated with CPU-time results. Application of Richardson extrapolation in time is then investigated to improve the first-order temporal accuracy, and finally, it is shown that numerical predictions agree with available experimental data during sub-picosecond laser heating of a thin film.