We develop a time-dependent theoretical numerical model to simulate the density profiles of the ions (i.e., O+(2P), O+(2D), N2 +, O+(4S), N+, O2 +, and NO+) and free electrons in E and F regions. In this model, the ion photoionization production rates, the photoelectron ionization production effect, and the chemical reactions between the ionized species and the neutral compositions are considered, and the plasma transport processes are not included. The simulation results show that the mean electron density ratios of the model simulations to the AE-C satellite measurements for the Solar Dynamics Observatory-Extreme Ultraviolet Variability Experiment (SDO-EVE), EUV flux model for Aeronomic Calculations, and Hinteregger-Fukui-Gilson solar irradiance models are, respectively, 0.97, 0.79, and 0.71 in a height range 140–400 km and 0.79, 0.75, and 0.64 in a height range 90–150 km. A comparison shows that the electron densities simulated by the model developed in this study are much more consistent with the in situ measurements made by the AE-C satellite than those predicted by the International Reference Ionosphere model and simulated by the Thermosphere Ionosphere Electrodynamics General Circulation Model. The model simulations with SDO-EVE solar irradiance input indicate that the photoelectron impact production process contributes about 20%–30% of the total atomic ion densities throughout the height range 130–400 km around noon. However, the photoelectron production effect on the molecular ion densities is very minor (less than about 7%) above 275 km. Below 250 km, its effect increases with the decrease of height for O2 + and N2 +, from about 4% and 10% at 250 km to 13% and 27% at 150 km, respectively.