The equation modelling the evolution of a foam (a complex porous medium consisting of a set of gas bubbles surrounded by liquid films) is solved numerically. This model is described by the reaction–diffusion differential equation with a free boundary. Two numerical methods, namely the fixed-point and the averaging in time and forward differences in space (the Crank–Nicolson scheme), both in combination with Newton’s method, are proposed for solving the governing equations. The solution of Burgers’ equation is considered as a special case. We present the Crank–Nicolson scheme combined with Newton’s method for the reaction–diffusion differential equation appearing in a foam breaking phenomenon.