Mathematical models of multi-step enzyme systems immobilized on porous electrodes are discussed. The model contains a non-linear term related to Monod kinetics. This work presents a better approximate analytical solution for non-linear differential equations of concentrations profile is derived using new approach of Homotopy perturbation method. Moreover, the derived analytical expressions are compared with numerical simulations. In addition, analytical expression for current density is derived and also it is compared with experimental data of already published work of Hettige et al [16]. The analytical expressions derived here are good enough to predict the behavior of the system. The influence of dimensionless governing parameters are discussed and presented graphically. The presented solutions are more reliable and easy to predict the dynamic behavior of the system by varying the parameters. © 2019, Springer Nature Switzerland AG.