Question: system of complicated ODEs

Dear guys! I want to solve a system of ODEs but I have some problems. Plaese have a look and let me know:

> Omega_m0 := 0.22: Zeta_m := 10^(-4): n := 0.35: h := 0.63: omega := 1:

> eq := z-> 1 = H^(2*n-2)*(1-Omega_m0)/(h^(2*n-2)) + h^2*Omega_m0*(1+z)^3*(1+Zeta_m*exp(-2*psi(z)))/(H^2) + omega*H^2*(1+z)^2*(diff(psi(z),z))^2/(6*h^2);

> yp := implicitdiff(eq(z), H, z);

> ode1 := diff(H(z),z)=subs(H=H(z),yp);

> ode2 := H(z)*(1+z)^2*diff(H(z),z)*diff(psi(z),z) - 2*H(z)^2*(1+z)*diff(psi(z),z) + H(z)^2*(1+z)^2*diff(diff(H(z),z),z) +6*Zeta_m*(1+z)^3*Omega_m0*h^2*exp(-2*psi(z))/omega = 0;

> sys := {ode1, ode2}:

> ics := {H(0) = 1, psi(0) = 0}:

> sol := dsolve(`union`(sys, ics), numeric, output = listprocedure, stiff = true);

The error I recieve is:  Error, (in dsolve/numeric/SC/IVPsetup) initial conditions must be numeric

I think my initial condition H(0)=1, is wrong. I think maybe this condition does not satisfy the system. If I am right how can find the right H(0)?

Thanks a lot.