Would it be possible, perhaps just to make your problem simple enough for me to understand, for you to expand your y'[t]=pflow[y[t]] into a list of explicit differential equations in terms of your y1[t], y2[t],...y6[t] and the derivatives of those and the coefficient variables that you have constructed? If that could be done then it might be possible for someone to see how to then translate that into a notation that NDSolve would correctly accept.

I think there might be a misunderstanding or miscommunication between what NDSolve is expecting and what pflow is providing.

If you look at the documentation for NDSolve you can see that it is expecting some sort of differential equation in the usual form and it will use one of the classical methods of numerically solving that.

It isn't clear to me how your list of n pairs of real numbers is supposed to supply that needed information. Can you perhaps explain what the meaning of the list of pairs is?

Reading your code, it seems that you might expect the list to provide y[t]. If that is the case then I am puzzled by your code.

Or perhaps you expect the list to provide the derivative of y[t] and you are wanting NDSolve to find y[t].

Suppose you could interpolate your data points and get a function that mostly goes through those points. Would that approximate y[t] or the derivative of y[t]? If it is the derivative then would the integral of that interpolation perhaps be what you are looking for?

If you can help readers get up to speed then perhaps someone can help. Thanks