Abstract A numerical procedure for the calculation of potential distributions in non-linear polar dielectrics, in spherical coordinates and with dirichlet boundary conditions, is presented. The basis of this iterative procedure is the difference equation obtained from the differential equation for the potential in non-linear polar dielectrics. The results obtained with the relaxation procedure outlined, for some particular cases, are shown, and also the comparison with the analytical results previously introduced by Coffey and Scaife, is presented. The numerical results agree very well with the analytical values, showing the correctness of the proposed numerical methods. It is shown that the numerical solution is more powerful than the analytical solution.