A hybrid finite difference method for the longitudinal dispersion equation, which is based on combining the Holly-Preissmann scheme with fifth-degree Hermite interpolating polynomial and the generalized Crank-Nicholson scheme, is described and comparatively evaluated with other characteristics-based numerical methods. Longitudinal dispersion of an instantaneously-loaded pollutant source is simulated, and computational results are compared with the exact solution. The present method is free from wiggles regardless of the Courant number, and exactly reproduces the location of the peak concentration. Overall accuracy of the computation increases for smaller value of the weighting factor, $ heta$of the model. Larger values of $ heta$ overestimates the peak concentration. Smaller Courant number yields better accuracy, in general, but the sensitivity is very low, especially when the value of $ heta$ is small. From comparisons with the hybrid method using cubic interpolating polynomial and with splitoperator methods, the present method shows the best performance in reproducing the exact solution as the advection becomes more dominant.