Dwarf spheroidal (dSph) galaxies are among the most promising targets for the indirect detection of dark matter (DM) from annihilation and/or decay products. Empirical estimates of their DM content - and hence the magnitudes of expected signals - rely on inferences from stellar-kinematic data. However, various kinematic analyses can give different results and it is not obvious which are most reliable. Using extensive sets of mock data of various sizes (mimicking ultra-faint and classical dSphs) and an MCMC engine, here we investigate biases, uncertainties, and limitations of analyses based on parametric solutions to the spherical Jeans equation. For a variety of functional forms for the tracer and DM density profiles, as well as the orbital anisotropy profile, we examine reliability of estimates for the astrophysical J- and D-factors for annihilation and decay, respectively. For large (N > 1000) stellar-kinematic samples typical of classical dSphs, errors tend to be dominated by systematics, which can be reduced through the use of sufficiently general and flexible functional forms. For small (N < 100) samples typical of ultrafaints, statistical uncertainties tend to dominate systematic errors and flexible models are less necessary. We define an optimal strategy that would mitigate sensitivity to priors and other aspects of analyses based on the spherical Jeans equation. We also find that the assumption of spherical symmetry can bias estimates of J (with the 95% credibility intervals not encompassing the true J-factor) when the object is mildly triaxial (axis ratios b/a = 0.8, c/a = 0.6). A concluding table summarises the typical error budget and biases for the different sample sizes considered.