In order to study beam optics of NIRS‐ECR ion source used in the HIMAC project, the EGUN code has been modified to make it capable of modeling ion extraction from a plasma. Two versions of the modified code are worked out with two different methods in which 1D and 2D sheath theories are used, respectively. Convergence problem of the strong nonlinear self‐consistent equations is investigated. Simulations on NIRS‐ECR ion source and HYPER‐ECR ion source are presented in this paper, exhibiting an agreement with the experiment results.