# ParametricNDSolve with more functions

A MWE of my problem is this:

Clear[k]
Ssol = ParametricNDSolveValue[{D[h00[t], {t, 2}] == -k h22[t],
D[h22[t], {t, 2}] == -k h00[t], h00[0.4] == 0, h00[0.99999] == 0,
h22[0.4] == 0, h22[0.99999] == 0}, {h00, h22}, {t, 0.4, 0.99999},
k]


This code evaluates (it seems) producing the following:

On the other hand the following code:

k = 1;
{S1, S2} =
NDSolveValue[{D[h00[t], {t, 2}] == -k h22[t],
D[h22[t], {t, 2}] == -k h00[t], h00[0.4] == 0, h00[0.99999] == 0,
h22[0.4] == 0, h22[0.99999] == 0}, {h00, h22}, {t, 0.4, 0.99999}]
Plot[{S1[t], S2[t]}, {t, 0.4, 0.9}]


produces a result with two functions which I can plot:

I expect that ParametricNDSolveValue to produce two interpolating functions, but it produces only one. It seems that ParametricNDSolveValue sees the two functions h00,h22 as one: how can I access to the solution for both functions? Or how can I fix my code?

If needed I use M. 11.3

ParametricNDSolve evaluates two interpolation-objects depending on the parameter k.
Ssol[1]

Through[Ssol[1][t]]