Nonlinear partial differential equations of the type
\[ u_t+\text{div}({\mathbf q}f(u)) -\Delta \Phi(u)=0 \]
are investigated especially from the numerical point of view. The function \(\Phi\) is supposed to be strictly increasing, but there can exist some values \(s\), such that \(\Phi^\prime(s)=0\). The numerical method is based on finite volumes and a new method for computing the numerical flux is proposed. The finite volume method is very convenient for a numerical approximation because these methods are suitable for purely nonlinear hyperbolic equations and they are efficient in the case of degenerate parabolic equations. The proposed method is a generalization of an idea of A. M. Il’in [Mat. Zametki 6, 237–248 (1969; Zbl 0185.42203)], D. N. de G. Allen and R. V. Southwell [Q. J. Mech. Appl. Math. 8, 129–145 (1955; Zbl 0064.19802)], D. L. Scharfetter and H. K. Gummel [Large signal analysis of a solicon Read diode. IEEE Transactions of Electron Devices 16, 64–77 (1969)] for the linear case in \(\Phi\) and \(f\). This method is based on setting the numerical flux \(G\) as a function \(G=g(a,b,q,h)\) where the function \(g\) is defined from the solution of the two point boundary value problem on \((0,h)\): \[ \left(-\Phi(v)^\prime+qf(v)\right)^\prime =0 \text{ on } (0,h) \]
\[ v(0)=a, \;\;v(h)=b. \] and \(g\) is set to be equal to the constant value \(\left(-\Phi(v(x))^\prime+qf(v(x))\right)^\prime\) for all \(x\in (0,h)\). The existence of a solution of such a Dirichlet two point boundary problem is proved. Then new numerical method including the numerical flux obtained by using the two point boundary method is proposed. Stability properties and finally convergence of the proposed scheme are proved. Efficiency and accuracy of the proposed scheme are confirmed by a practical numerical experiment where a comparison of the proposed method and other numerical schemes is presented.