Finding the magnetic flux mapping in the ionosphere is very important. It would not only divide the surface into the elements with the same flux, but also indicate locations of conjugated points. It is important for studies of field aligned currents and bouncing of energetic charged particles and their precipitation. The existing methods involve numerical magnetic field lines tracing in the entire volume of the magnetosphere or numerical integration along assumed contour lines of the Euler potentials on the surface of the ionosphere. It is possible to determine the mapping with these methods near the magnetic equator, but not on middle latitudes and near and inside the polar caps. Our approach is to search for the Euler potentials as a sum of basic functions with their coefficients. Each basic function is a product of a sine or cosine of longitude multiplied by m and the Legendre polynomial of the colatitude angle cosine and of the order n. Maxima of m and n in this calculation were set to 13. The difference between the radial component from the cross product of the Euler potentials gradients and from International Geomagnetic Reference Field is less than 0.01 percent. We discuss the possibility of using orthogonal coordinates defined on the sphere's surface, which remain finite functions of θ and φ everywhere except for the vicinities of the North and South poles. The issues with numerical errors accumulated on long tracing are avoided when using this approach.