sorry about the formatting. I don't know why it cut off the first part of the expressions.

It should read (just in ascii now)

dsys:={D(a)(t)=sqrt(RhoDE(phi(t),a(t))+ RhoDM(phi(t),a(t)) + RhoR(phi(t),a(t)))*a(t),

D(D(phi))(t) + 3*a(t)*D(phi)(t)/D(a)(t)*dtOm(phi(t),a(t))=0,

a(1000)=#

phi(1000)=#,

D(phi)(1000)=#,

}

(the ICs are given as expressions  in my worksheet and above but they evaluate numerically, so I just put # this time)

Where I've already defined the functions

dtOM:= (phi,a)--> -alpha*M^(alpha+4)/phi^(alpha+1) + 1/(3*pi^2)*integrate(2*g^2*phi^2*k^4*(a*sqrt*(g^2*phi^2+k^2)/T+1)*exp(-a*sqrt*(g^2*phi^2+k^2)/T)*(g^2*phi^2+k^2)^(-3/2), k, 0..infinity)

where M, g, alpha and T are constants.

and

RhoDE:= (phi,a)-->M^(4(alpha+1))/phi^alpha+D(phi)(t)/2*a^2

RhoDM:=a--> C/a^3

RhoR:=pi^2/15*(T/a)^4

with C another constant.

Like I sad above, when I call dsolve, the error I get is a bit confusing.

>dsolve(dsys, numeric)

Error, (in unknown) invalid input: `dsolve/numeric/lexorder` expects its 2nd argument, b, to be of type {indexed, name, numeric}, but received phi(t)

Any help?

Thanks for reading after that disasterous first post. Sorry about that.