A simple Mathematica (version 7) code for computing S-state energies and wave functions of two-electron (helium-like) ions is presented. The elegant technique derived from the classical papers of Pekeris (1958, 1959, 1962, 1965, 1971) [1–3] is applied. The basis functions are composed of the Laguerre functions. The method is based on the perimetric coordinates and specific properties of the Laguerre polynomials. Direct solution of the generalized eigenvalues and eigenvectors problem is used, distinct from the Pekeris works. No special subroutines were used, only built-in objects supported by Mathematica. The accuracy of the results and computation times depend on the basis size. The ground state and the lowest triplet state energies can be computed with a precision of 12 and 14 significant figures, respectively. The accuracy of the higher excited states calculations is slightly worse. The resultant wave functions have a simple analytical form, that enables calculation of expectation values for arbitrary physical operators without any difficulties. Only three natural parameters are required in the input. The above Mathematica code is simpler than the earlier version (Liverts and Barnea, 2010 [4]). At the same time, it is faster and more accurate. Program summary Program title: TwoElAtomSL(SH) Catalogue identifier: AEHY_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEHY_v1_0.html Program obtainable from: CPC Program Library, Queenʼs University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 11 434 No. of bytes in distributed program, including test data, etc.: 540 063 Distribution format: tar.gz Programming language: Mathematica 7.0 Computer: Any PC Operating system: Any which supports Mathematica; tested under Microsoft Windows XP and Linux SUSE 11.0 RAM: ⩾ 10 9 bytes Classification: 2.1, 2.2, 2.7, 2.9 Nature of problem: The Schrödinger equation for atoms (ions) with more than one electron has not been solved analytically. Approximate methods must be applied in order to obtain the wave functions or another physical attributes from quantum mechanical calculations. Solution method: The S-wave function is expanded into a triple set of basis functions which are composed of the exponentials combined with the Laguerre polynomials in the perimetric coordinates. Using specific properties of the Laguerre polynomials, solution of the two-electron Schrödinger equation reduces to solving the generalized eigenvalues and eigenvector problem for the proper Hamiltonian. The unknown exponential parameter is determined by means of minimization of the corresponding eigenvalue (energy). Restrictions: First, the too large length of expansion (basis size) takes the too large computation time and operative memory giving no perceptible improvement in accuracy. Second, the number of shells Ω in the wave function expansion enables one to calculate the excited nS-states up to n = Ω + 1 inclusive. Running time: 2–60 minutes (depends on basis size and computer speed).