We propose a numerical method for the simulation of a quasi-linear parabolic biofilm model that exhibits three non-linear diffusion effects: (i) a power law degeneracy, (ii) a super diffusion singularity and (iii) non-linear cross-diffusion. The method is based on a spatial Finite Volume discretisation in which cross-diffusion terms are formally treated as convection terms. Time-integration of the resulting semi-discretised system is carried out using an error-controlled, time-adaptive, embedded Rosenbrock–Wanner method. We compare several variants of the method and two variants of the model to investigate how details such as the choice cross-diffusion coefficients, and specific variants of the time integrator affect simulation time.