Question: Numerical solution of ode system

Hi,

I have been trying to get a numerical solution to the following ode system for weeks. I am new to Maple, and I think my issues might be with the coding rather than the maths.

I get the error message:

Error, (in dsolve/numeric/checksing) ode system has a removable singularity at t=0. Initial data is restricted to {v(t) = 0., Phi(t) = 0.}

I know there is a singularity at t=0, so I tried setting the 'range' option to something just above 0, but it is still not working.

Here is my worksheet (without the output, and I am using Maple 10 classic). Sorry, it is pretty ugly:

restart;with(DEtools):with(plots):with(linalg):infolevel[dsolve]:=4:Digits:=20;
> k:=0.065;
> h:=0.5;
> H0:=100*h;
> K:=(omm,omv)->-H0^2*(1-omm-omv);
> Lambda:=(omv)->3*H0^2*omv;
> G:=6.673e-11;#m^3 kg^-1 sec^-2
> a0:=(omm)->25000*omm*h^2:
> a0(1);
> adotarad:=(t)->1/t;
> adotamat:=(t)->(2/t);
t = conformal time:
> arad:=(t)->t/a0(.3):
> amat:=(t)->t^(2)/a0(.3):
> raddom:=loglogplot(arad(t),t=1e-6..1,axes=framed):
> matdom:=loglogplot(amat(t),t=1..a0(.3)^(3/2),axes=framed):
> #display(raddom,matdom);
> rho_dm:=0.3: rho_r:=0.7:

> de1:=diff(Theta0(t),t)=-k*Theta1(t)+(3*adotamat(t))^(-1)*(k^2+3*adotamat(t)^(2)*Phi(t)-4*evalf(Pi)*G*(rho_dm*delta(t)+ 4*rho_r*Theta0(t))):
> de2:=diff(Theta1(t),t)=-k/3*Phi(t)+k/3*Theta0(t):
> de3:=diff(delta(t),t)=-adotamat(t)^(-1)*(k^2+3*adotamat(t))^2*Phi(t)-4*evalf(Pi)*G*rho_dm*delta(t)+ 16*evalf(Pi)*G*rho_dm**rho_r*Theta0(t)-I*k*v(t):
> de4:=diff(v(t),t)=I*k*Phi(t)-adotamat(t)*v(t):
>de5:=diff(Phi(t),t)=3*adotamat(t)*(4*evalf(Pi)*G*amat(t)^2*rho_dm*delta(t)+16*evalf(Pi)*G*amat(t)^2*rho_r*Theta0(t))-(k^2+3*adotamat(t)^2)*Phi(t):

> deqs := de1,de2,de3,de4,de5;
> ics := Theta0(0)=1e-3;Theta1(0)=1e-3;Phi(0)=0.002;v(0)=0.1;delta(0)=0.003;
> #ics := Phi(0)=0.1;Theta0(0)=0.001,Theta1(0)=0.001,delta(0)=0.001,v(0)=0.1;
> fns := Phi(t),Theta0(t),Theta1(t),delta(t),v(t);
> eqsol := dsolve({deqs,ics},{fns},numeric,'range'=1..1e6);

Any help/comments hugely appreciated.

Peter