A study of the $B_c^+\to J/\psi D_s^+$ and $B_c^+\to J/\psi D_s^{*+}$ decays using 139\,fb$^{-1}$ of integrated luminosity collected with the ATLAS detector from $\sqrt{s} = 13$\,TeV $pp$ collisions at the LHC is presented. The ratios of the branching fractions of the two decays to the branching fraction of $B_c^+\to J/\psi \pi^+$ decay channel are measured: $\mathcal B(B_c^+\to J/\psi D_s^+)/\mathcal B(B_c^+\to J/\psi \pi^+) = 2.76\pm 0.33\pm 0.29\pm 0.16$, $\mathcal B(B_c^+\to J/\psi D_s^{*+})/\mathcal B(B_c^+\to J/\psi \pi^+) = 5.33\pm0.61\pm 0.67\pm 0.32$. Here the first error is statistical, the second is systematic, and the third is due to uncertainty on the $D_s^+\to\phi(K^+ K^-)\pi^+$ decay branching fraction. The ratio between the branching fractions of the two decays is found to be $\mathcal B(B_c^+\to J/\psi D_s^{*+})/\mathcal B(B_c^+\to J/\psi D_s^+) = 1.93\pm0.24 \pm 0.09$. For the $B_c^+\to J/\psi D_s^{*+}$ decay, the transverse polarisation fraction, $\Gamma_{\pm\pm}/\Gamma$, is measured to be $0.70\pm0.10\pm 0.04$. The precision of the measurements exceeds that of all previous studies of these decays.