Axisymmetric equilibrium of a magnetically confined plasma with multiple ion species is considered. To describe hot plasmas with isothermal surfaces, we adopt a formulation consistent with zero poloidal ion flow. This formulation includes all magnetic field components and also correctly includes all effects of toroidal ion rotation. There are two free surface functions for each species and a third which is determined by a differential equation relating surface functions per species. We have developed and implemented an algorithm for the solution of the resulting nonlinear equations and found solutions with large charge and mass contrast among the ion species for both compact (r = 0 included) and annular (r = 0 excluded) domains. Our solution method allows for arbitrary domain shapes, includes far-field conditions, and treats any combination of electrically conducting or insulating walls. Appropriate surface functions are used to describe the transition from closed to open field plasma in a reasonable manner. Solutions for advanced fuel cycle fusion systems (both D–3He and p–11B) are presented to illustrate the power of the method. Finally, we briefly discuss the special issues associated with obtaining very elongated solutions and describe the algorithm for implementing these features.