In this paper, a nonuniform size modified Poisson-Boltzmann equation (SMPBE) for a protein in a solvent with multiple ionic species in distinct ionic sizes is derived by using a new electrostatic free energy functional and solution decomposition techniques. It is then proved to have a unique solution, and the solution satisfies a system consisting of nonlinear algebraic equations and one Poisson dielectric interface problem. To solve it numerically, two new finite element iterative schemes are proposed by using nonlinear successive over-relaxation techniques, along with an improved uniform SMPBE for generating initial iterates. Furthermore, they are programmed in Python and Fortran as a software package for solving the nonuniform SMPBE, and numerically tested on a Born ball test model and a protein in a sodium chloride solution and a sodium chloride and potassium chloride solution. Numerical results confirm the convergence of the two new iterative schemes and demonstrate the high performance of the new software package.