Using ParametricNDSolveValue with initial conditions of another ParametricNDSolveValue

Question

I have a system of parameter-dependent ODEs, which I can solve using ParametricNDSolveValue. Now, I have another ODE, which has as initial condition the solution of one component of the first system at a specific time point.

Consider the example

pfun = ParametricNDSolveValue[{x'[t] == a*x[t], x[0] == 1}, {x}, {t, 0, 10}, {a}];
qfun = ParametricNDSolveValue[{y'[s] == a^2*y[s], y[0] == First[pfun[a]][2]}, {y}, {s, 0, 10}, {a}];

Then,

qfun[4]

gives the output: ParametricNDSolveValue::ndinnt: "Initial condition 4.[2.] is not a number or a rectangular array of numbers".

This is odd, since First[pfun[.5]][2] produces $2.71828$ as output, as one expects.

However, omitting the curly brackets of x in the first line gives the correct result:

pfun = ParametricNDSolveValue[{x'[t] == a*x[t], x[0] == 1}, x, {t, 0, 10}, {a}];
qfun = ParametricNDSolveValue[{y'[s] == a^2*y[s], y[0] == pfun[a][2]}, {y}, {s, 0, 10}, {a}];
First[qfun[4]][2]

produces the output $2.35385*10^{17}$, which is $\exp(2a)\exp(a^2s)$ for $a=4$ and $s=2$.

As mentioned, pfun is a system of equations in the real problem, so omitting curly brackets is no option for me.

Any ideas?

UPDATE

I found out that the problem with the first two lines of code is that the fragment First[pfun[a]] is first evaluated with a having no value, thus First[pfun[a]] reduces to a. How can I specify that First is evaluated after a is given a specific value?

Answer 1

this looks a little like the more common case where you want to evaluate only for numeric arguments. Here you want to only evaluate for lists, here is one way to achieve that:

pfun = ParametricNDSolveValue[{x'[t] == a*x[t], x[0] == 1}, {x}, {t,0, 10}, {a}];

first[l_List] := First[l];

qfun = ParametricNDSolveValue[
   {y'[s] == a^2*y[s],y[0] == first[pfun[a]][2]}, {y}, {s, 0, 10}, {a}
];

qfun[4]

for your more general case of a system of ODEs you might want to define something like:

part[l_List, n_Integer] := Part[l, n];

and then use it like:

qfun = ParametricNDSolveValue[
   {y'[s] == a^2*y[s], y[0] == part[pfun[a], 1][2]}, {y}, {s, 0, 10}, {a}
];