How to integrate ParametricNDSolve solution with respect to a parameter?

Question

I have just upgraded to the new version of Mathematica because of its new built-in ParametricNDSolve function. I need to solve a first-order non-linear ordinary differential equation that has two parameters. I have no problem obtaining appropriate solutions to this equation, but once I have the solution I need to integrate out the dependence on the parameters. For instance, I would like to write something like

sol=ParametricNDSolveValue[{-y*a'[x]==(a[x]^2-1)f[x]+a[x],a[-10]==0.6},a,{x,-10,10},{y}]
Integrate[sol[y][x],{y,0,1}]

If I have an InterpolationFunction which has two arguments I can very easily integrate over one of the two arguments, but this does not seem possible to do with a ParametricFunction. The best idea I have had so far is to form an InterpolationFunction from the ParametricFunction and then integrate from there, but that seems inefficient and I would like something a little faster and more straightforward. Any help would be greatly appreciated.

Answer 1

Here is a way:

sol = ParametricNDSolveValue[{-y*a'[x] == (a[x]^2 - 1)*x + a[x], 
   a[-10] == 0.6}, a, {x, -10, 10}, {y}]

sol[1]

NIntegrate[sol[y][t] /. y -> 1, {t, 0, 1}]

or alternatively via a Table:

ListLinePlot[
 Table[NIntegrate[sol[y][t], {t, 0, 1}], {y, 0.1, 1, 0.1}]]