We present new dynamical models of dwarf spheroidal galaxies (dSphs) in which both the stellar component and the dark halo are described by analytic distribution functions that depend on the action integrals. In their most general form these distribution functions can represent axisymmetric and possibly rotating stellar systems. Here, as a first application, we model the Fornax dSph, limiting ourselves, for simplicity, to the non rotating, spherical case. The models are compared with state-of-the-art spectroscopic and photometric observations of Fornax, exploiting the knowledge of the line-of-sight velocity distribution of the models and accounting for the foreground contamination from the Milky Way. The model that best fits the structural and kinematic properties of Fornax has a cored dark halo, with core size $r_{rm c}simeq1.03$ kpc. The dark-to-luminous mass ratio is $(M_{rm dm}/M_{star})|_{R_{rm eff}}simeq9.6$ within the effective radius $R_{rm eff} simeq 0.62,$kpc and $(M_{rm dm}/M_{star})|_{3 {rm kpc}} simeq 144$ within 3 kpc. The stellar velocity distribution is isotropic almost over the full radial range covered by the spectroscopic data and slightly radially anisotropic in the outskirts of the stellar distribution. The dark-matter annihilation $J$-factor and decay $D$-factor are, respectively, $log_{10}(J$ $[$GeV$^2$ cm$^{-5}])simeq18.34$ and $log_{10}(D$ $[$GeV cm$^{-2}])simeq18.55$, for integration angle $theta = 0.5^{circ}$. This cored halo model of Fornax is preferred, with high statistical significance, to both models with a Navarro, Frenk and White dark halo and simple mass-follows-light models.