A biofilm model is formulated by a quasilinear reaction–diffusion system of equations. It exhibits possible degeneracy and singularities in the diffusion coefficient. An original combination of Crank–Nicolson and Newton methods speed up computations. Numerical tests show the expected behaviour of biomass and nutrients concentrations. Numerical results also illustrate some theoretically proved qualitative properties.