This paper is devoted to the study of the large time behaviour of viscosity solutions of parabolic equations with Neumann boundary conditions. This work is the sequel of Hu et al. (SIAM J Control Optim 53:378-398, 2015) in which a probabilistic method was developed to show that the solution of a parabolic semilinear PDE behaves like a linear term shifted with a function v, where is the solution of the ergodic PDE associated to the parabolic PDE. We adapt this method in finite dimension by a penalization method in order to be able to apply an important basic coupling estimate result and with the help of a regularization procedure in order to avoid the lack of regularity of the coefficients in finite dimension. The advantage of our method is that it gives an explicit rate of convergence.