  We describe numerical tools for the stability analysis of extrasolar
planetary systems. In particular, we consider the relative Poincare variables
and symplectic integration of the equations of motion. We apply the tangent map
to derive a numerically efficient algorithm of the fast indicator MEGNO (a
measure of the maximal Lyapunov exponent) that helps to distinguish chaotic and
regular configurations. The results concerning the three-planet extrasolar
system HD 37124 are presented and discussed. The best fit solutions found in
earlier works are studied more closely. The system involves Jovian planets with
similar masses. The orbits have moderate eccentricities, nevertheless the best
fit solutions are found in dynamically active region of the phase space. The
long term stability of the system is determined by a net of low-order two-body
and three-body mean motion resonances. In particular, the three-body resonances
may induce strong chaos that leads to self-destruction of the system after Myrs
of apparently stable and bounded evolution. In such a case, numerically
efficient dynamical maps are useful to resolve the fine structure of the phase
space and to identify the sources of unstable behavior.
