I am trying to use finite difference method to solve for $u(x,t)$ in the equation:

\begin{align} \frac{\partial^2u}{\partial x^2} = \frac{au}{1+bu}, \end{align} which is actually part of a system of PDEs. The equation came from the Michaelis-Menten law used in modelling tumor growth where $u(x,t)$ is the oxygen tension. The RHS of the equation can be re-written as follows, \begin{align} \frac{au}{1+bu} = \frac{a}{b}\left(1-\frac{1}{1+bu}\right). \end{align} The boundary conditions are $\partial u/\partial x = 0$ at $x=0$ and $u(1,t)=0$. Usually, for Poisson equation $\dfrac{\partial^2u}{\partial x^2}=f(x)$, which is quite similar to the above, I just do \begin{align} \frac{u_{i+1}-2u_i+u_{i-1}}{\Delta x^2}=f_i, \end{align} for each $i=1,...,n$ and re-write the resulting discretised equations in a matrix form \begin{align} {\bf Au=b}, \end{align} and thus manipulate the matrix system (or use any other methods like Gauss-Seidel) to get the solution ${\bf u}$. But, how do I do it for the equation above? I would be very grateful if someone could give me a small hint or point me to some useful reference.