We develop numerical methods for solving nonlinear equations of conservation laws with flux function that depends on discontinuous coefficients. Using a relaxation approximation, the nonlinear equation is transformed to a semilinear diagonalizable problem with linear characteristic variables. Eulerian and Lagrangian methods are used for the advection stage while an implicit