# How to integrate ParametricNDSolve solution with respect to a parameter?

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.

-
it's always a good idea to post code that works when copying it into a fresh session. – user21 Apr 29 '13 at 6:21

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}]]


-