In a recently published article of Daripa and Pasa [Transp. Porous Media (2007) 70:11-23], the stabilizing effect of diffusion in three-layer Hele-Shaw flows was proved using an exact analysis of normal modes. In particular, this was established from an upper bound on the growth rate of instabilities which was derived from analyzing stability equations. However, the method used there is not constructive in the sense that the upper bound derived from actual numerical discretization of the problem could be significantly different from the exact one reported depending on the scheme used. In this paper, a numerical approach to solve the stability equations using a finite difference scheme is presented and analyzed. An upper bound on the growth rate is derived from numerical analysis of the discrete system which also shows the diffusive slowdown of instabilities. Upper bounds obtained by this numerical approach and by the analytical approach are compared. The present approach is constructive and directly leads to the implementation of the numerical approach to obtain approximate solutions in the presence of diffusion. The contributions of the paper are the novelty of the approach and a bound on the growth rates that does not depend on the solution itself.