By means of polynomial approximation and iteration procedure using linearization solution as initial solution, the non-linear Poisson–Boltzmann equation describing a spherical colloidal particle immersed in an arbitrary valence and mixed electrolyte solution is solved analytically, and analytical expressions for electrical potential distribution ψ(r) and surface charge density/surface potential relationship (σ/ψ 0) are acquired. The σ/ψ 0 expression performs very well for the entire range of a reduced colloidal radius x 0 only if a reduced surface potential y 0 is lower than 15.5; particularly, in the case of x 0 >10, the σ/ψ 0 expression applies over the entire range of y 0, and maximum of the absolute value of percent relative error (PRE) is not larger than 2 if one does not come across a parameter combination as of an extreme case of x 0 < 0.9 and y 0 > 15.5. The ψ(r) expression outperforms the linearization solution greatly and performs very well for the domain of x 0 ≤ 1.5 and y 0 ≤ 4.5 with maximum of the absolute value of PRE not larger than 3, whereas that of the linearization solution may stand 22 for the same parameter range. It is concluded that the present σ/ψ 0 expression is the first one valid for entire range of x 0, given that y 0 is not unreasonably high, and the present paper proposes, to some extent, a “universal” way for solving nonlinear differential equations.