A steady state Poisson-Nernst-Planck (PNP) system is studied both analytically and numerically with particular attention on I-V relations of ion channels. Assuming the dielectric constant $\varepsilon$ is small, the PNP system can be viewed as a singularly perturbed system. Due to the special structures of the zeroth order inner and outer systems, one is able to derive more explicit expressions of higher order terms in asymptotic expansions. For the case of zero permanent charge, under the assumption of electro-neutrality at both ends of the channel, our result concerning the I-V relation for two oppositely charged ion species is that the third order correction is \textit{cubic} in $V$, and, furthermore (Theorem \ref{3rd order}), up to the third order, the cubic I-V relation has \textit{three distinct real roots} (except for a very degenerate case) which corresponds to the bi-stable structure in the FitzHugh-Nagumo simplification of the Hodgkin-Huxley model. Three numerical experiments are conducted to check the cubic-like feature of the I-V curve, study the boundary value effect on the I-V relation and investigate the permanent charge effect on the I-V curve, respectively.