The steady-state Poisson-Nernst-Planck (ssPNP) equations are an effective model for the description of ionic transport in ion channels. It is observed that an ion channel exhibits voltage-dependent switching between open and closed states. Different conductance states of a channel imply that the ssPNP equations probably have multiple solutions with different level of currents. We propose numerical approaches to study multiple solutions to the ssPNP equations with multiple ionic species. To find complete current-voltage (I-V) and current-concentration (I-C) curves, we reformulate the ssPNP equations into four different boundary value problems (BVPs). Numerical continuation approaches are developed to provide good initial guesses for iteratively solving algebraic equations resulting from discretization. Numerical continuations on V, I, and boundary concentrations result in S-shaped and double S-shaped (I-V and I-C) curves for the ssPNP equations with multiple species of ions. There are five solutions to the ssPNP equations with five ionic species, when an applied voltage is given in certain intervals. Remarkably, the current through ion channels responds hysteretically to varying applied voltages and boundary concentrations, showing a memory effect. In addition, we propose a useful computational approach to locate turning points of an I-V curve. With obtained locations, we are able to determine critical threshold values for hysteresis to occur and the interval for V in which the ssPNP equations have multiple solutions. Our numerical results indicate that the developed numerical approaches have a promising potential in studying hysteretic conductance states of ion channels.