That's the whole issue, that NDSolve, just like you, wants to see an explicit form of the RHS of the DE p'[t]==pflow[p], but I dont have one, since the function pflow[p] is given by some complicated program, which cannot be summarized by a "formula" (probably it can, in my case, with some effort, but it's very unnatural, and one should not be obliged to do it, the program should serve me, not the other way around). In particular, pflow can only accept as input vectors p with numerical entries, not abstract symbols.

OK, so I can tell you that after searching around and asking some experts, I learnt that there is a very simple solution to this problem in Mathematica, but which is not obvious at all and is not mentioned in the documentation (at least I couldnt find it). You just add the single word "List" after the argument in the definition of pflow (the RHS of the DE)

pflow[p_List] := Module[{...}, etc

That's all. After this you can write

psol = NDSolveValue[{p'[t] == pflow[p[t]], p[0] == p0}, p, {t, 0, T}];

and all goes well. Maybe it's in the documentation somewhere, but I didnt find it. It somehow tells NDSolve, when trying to acces pflow before starting the integration process, not to use it prematurely, but wait patiently, start with the available initial conditions which turn p into a List of numbers (without them p is just a symbol without any meaning), which will be then accepted by pflow and produce, using the integrating scheme used by NDSolve, the next value of p, etc.

Thanks again for trying to help.