  Advantageous numerical methods for solving the Dirac equations are derived.
They are based on different stochastic optimization techniques, namely the
Genetic algorithms, the Particle Swarm Optimization and the Simulated Annealing
method, their use of which is favored from intuitive, practical, and
theoretical arguments. Towards this end, we optimize appropriate parametric
expressions representing the radial Dirac wave functions by employing methods
that minimize multi parametric expressions in several physical applications. As
a concrete application, we calculate the small (bottom) and large (top)
components of the Dirac wave function for a bound muon orbiting around a very
heavy (complex) nuclear system (the $^{208}$Pb nucleus), but the new approach
may effectively be applied in other complex atomic, nuclear and molecular
systems.
