  By solving the Muskhelishvili-Omn\`es integral equations, the scalar form
factors of the semileptonic heavy meson decays $D\to\pi \bar \ell \nu_\ell$,
$D\to \bar{K} \bar \ell \nu_\ell$, $\bar{B}\to \pi \ell \bar\nu_\ell$ and
$\bar{B}_s\to K \ell \bar\nu_\ell$ are simultaneously studied. As input, we
employ unitarized heavy meson-Goldstone boson chiral coupled-channel amplitudes
for the energy regions not far from thresholds, while, at high energies,
adequate asymptotic conditions are imposed. The scalar form factors are
expressed in terms of Omn\`es matrices multiplied by vector polynomials, which
contain some undetermined dispersive subtraction constants. We make use of
heavy quark and chiral symmetries to constrain these constants, which are
fitted to lattice QCD results both in the charm and the bottom sectors, and in
this latter sector to the light-cone sum rule predictions close to $q^2=0$ as
well. We find a good simultaneous description of the scalar form factors for
the four semileptonic decay reactions. From this combined fit, and taking
advantage that scalar and vector form factors are equal at $q^2=0$, we obtain
$|V_{cd}|=0.244\pm 0.022$, $|V_{cs}|=0.945\pm 0.041$ and $|V_{ub}|=(4.3\pm
0.7)\times10^{-3}$ for the involved Cabibbo-Kobayashi-Maskawa (CKM) matrix
elements. In addition, we predict the following vector form factors at $q^2=0$:
$|f_+^{D\to\eta}(0)|=0.01\pm 0.05$, $|f_+^{D_s\to K}(0)|=0.50 \pm 0.08$,
$|f_+^{D_s\to\eta}(0)|=0.73\pm 0.03$ and $|f_+^{\bar{B}\to\eta}(0)|=0.82 \pm
0.08$, which might serve as alternatives to determine the CKM elements when
experimental measurements of the corresponding differential decay rates become
available. Finally, we predict the different form factors above the
$q^2-$regions accessible in the semileptonic decays, up to moderate energies
amenable to be described using the unitarized coupled-channel chiral approach.
