We develop a code to compute different chaos indicators for the planetary N body problem. The code uses a second order symplectic integrator to advance the Hamiltonian flow of the system in heliocentric canonical coordinates, and the corresponding variational equations. The chaos indicators are computed over a dense grid of initial conditions, allowing to map the stability of the phase space in detail. The code computes the MEGNO indicator and the time variance of the orbital elements, among others. We take advantage of the GPU parallelization with CUDA to simulate the nodes of the grid simultaneously. We apply this code to analyze the stability of the Kepler-46 system of exoplanets. This system is constituted by a Jupiter and a Saturn size planets in tightly packed orbits. The presence of a third, super-Earth size, planet has been also hypothesized. Our analysis helps to reveal the main resonant structure of the system, which shows many similarities with our own solar system. All the planets in the system lie in regions of regular motion, even considering their orbital uncertainties. Our results also put constraints to the possible radial migration of the planets in the system, indicating that the pair of Jovian planets could have never got closer than the location of the 2:3 resonance.