Traditional unsaturated flow models use a capillary pressure-saturation relationship determined under static conditions. Recently it was proposed to extend this relationship to include dynamic effects and in particular flow rates. In this paper, we consider numerical modeling of unsaturated flow models incorporating dynamic capillary pressure terms. The resulting model equations are of nonlinear degenerate pseudo-parabolic type with or without convection terms, and follow either Richards' equation or the full two-phase flow model. We systematically study the difficulties associated with numerical approximation of such equations using two classes of methods, a cell-centered finite difference method (FD) and a locally conservative Eulerian-Lagrangian method (LCELM) based on the finite difference method. We discuss convergence of the methods and extensions to heterogeneous porous media with different rock types. In convection-dominated cases and for large dynamic effects instabilities may arise for some of the methods while those are absent in other cases.