Thanks very much for your response Doug. The point about trying to start things off in a singularity is clear to me now!

Yes, the solutions should be complex-valued. They represent the evolution of matter density perturbations in the universe, and these oscillated in the early universe when things were much hotter.

I need to think about the time scale issue now.  It is a great help just to have the numerical integration working.

One more quick question: does odeplot have a logarithmic axes option?

Peter