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 [13] in which a probabilistic method was developped to show that the solution of a parabolic semilinear PDE behaves like a linear term $λT$ shifted with a function $v$, where $(v,λ)$ 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.