If secondary hydrocarbon recovery methods fail because of the occurrence of gravity override or viscous fingering one can turn to an enhanced oil recovery method like the injection of foam. The generation of foam can be described by a set of partial differential equations with strongly nonlinear functions, which impose challenges for the numerical modeling. To analyze the effect of foam on viscous fingering, we study the dynamics of a simple foam model based on the Buckley-Leverett equation. Whereas the Buckley-Leverett flux is a smooth function of water saturation, the foam will cause a rapid increase of the flux function over a very small saturation scale. Consequently its derivatives can become extremely large and impose a severe constraint on the time step due to the CFL condition. Until now, the methods applied to foam EOR processes are only first-order accurate and do not incorporate stabilization near the foam front as far as we know. In order to improve the accuracy near the foam front we make use of total variation diminishing schemes that preserve the numerical stability of the solution. Two dimensional simulations, including gravity, will shed light on the conditions under which foam might exhibit viscous fingering behavior.