Thank you very much Bill for trying to help.
The data structure that the function pflow uses, both for input and output, is a nested list of n pairs of real numbers (a vector of n points in R^2). For example, for n=3, something like {{1,2}, {3,4}, {5,6}}. I tested it by inserting the command pflow[p0] before the NDsolve command, and it seems to do its job properly. But somehow, when NDsolve calls pflow, something goes wrong.
Thanks again for yr effort, I really appreciate it.
