# numerical integration of equation obtained by ParametricNDSolve

I am trying to find a definite integral of the solution obtained using ParametricNDSolve (i.e. y[x,a,b]) from 0 to 1 and plot the obtained solution in two variables (i.e. a,b) with a. Can you please suggest me a solution for the same as my equation is given by

ParametricNDSolve[{Cos[y[x]]*y''[x]*y'[x] + Sin[y[x]]*y'''[x] -
Sin[y[x]]*y'[x]*y'''[x] + Cos[y[x]]*y''''[x] +
4*a*(Sin[y[x]])*(Cos[y[x]])^3 - a*(Sin[y[x]])*(Cos[y[x]])^2 +
b*Sin[y[x]]*Cos[y[x]] == 0, y[0] == 0, y''[1] == 0, y'[1] == 0,
y'''[1] == -1*a*(Cos[y[1]])}, y, {x, 0, 1}, {a, b}]


all I want is a plot of the integral of Sin[y[x]] from 0 to 1 with a and b.

One can integrate an expression within an NDSolve/ParametricNDSolve command. To integrate Sin[y[x]], use the DE u'[x] == Sin[y[x]].

psol = ParametricNDSolve[
{Cos[y[x]]*y''[x]*y'[x] + Sin[y[x]]*y'''[x] -
Sin[y[x]]*y'[x]*y'''[x] + Cos[y[x]]*y''''[x] +
4*a*(Sin[y[x]])*(Cos[y[x]])^3 - a*(Sin[y[x]])*(Cos[y[x]])^2 +
b*Sin[y[x]]*Cos[y[x]] == 0, y[0] == 0, y''[1] == 0, y'[1] == 0,
y'''[1] == -1*a*(Cos[y[1]]),
u'[x] == Sin[y[x]], u[0] == 0},
{y, u}, {x, 0, 1}, {a, b}];

Plot3D[u[a, b][1] /. psol // Evaluate, {a, 0, 3}, {b, 0, 3}]


I am not sure I understand your question correctly. Do these commands do what you want?

Plot[y[2, 3][x] /. sol, {x, 0, 1}]

plPoints =
Table[With[{f = Sin[y[a, b][x]] /. sol},
NIntegrate[f, {x, 0, 1}]], {a, 0, 3, 0.1}, {b, 0, 3, 0.1}]

ListPointPlot3D[plPoints]