The size-modified Poisson-Boltzmann equation (SMPBE) has been developed to consider ionic size effects in the calculation of electrostatic potential energy, but its numerical solution for a biomolecule remains a challenging research issue. To address this challenge, in this paper, we propose a solution decomposition formula and then develop an effective minimization protocol for solving the nonlinear SMPBE model by using finite element approximation techniques. As an application, a particular SMPBE numerical algorithm is constructed and programmed as a finite element program package for a biomolecule (e.g., protein and DNA) in a symmetric 1:1 ionic solvent. We also construct a nonlinear SMPBE ball model with an analytical solution and use it for validation of our new SMPBE numerical algorithm and program package. Furthermore, numerical experiments are made on a central charged ball model to show some physical features captured by the SMPBE model. Finally, they are made for six biomolecules with different net charges to demonstrate the computer performance of our SMPBE finite element program package. Numerical results show that the SMPBE model can capture some physical properties of an ionic solvent more reasonably, and can be solved more efficiently than the classic PBE model. As an application of the SMPBE model, free solvation energies were calculated and compared to the case of the PBE model.