  By applying the recurrence approach and coupling constant metamorphosis, we
construct higher order integrals of motion for the Stackel equivalents of the
$N$-dimensional superintegrable Kepler-Coulomb model with non-central terms and
the double singular oscillators of type ($n, N-n$). We show how the integrals
of motion generate higher rank cubic algebra $C(3)\oplus L_1\oplus L_2$ with
structure constants involving Casimir operators of the Lie algebras $L_1$ and
$L_2$. The realizations of the cubic algebras in terms of deformed oscillators
enable us to construct finite dimensional unitary representations and derive
the degenerate energy spectra of the corresponding superintegrable systems.
