Question: Having trouble solving my system of differential equations...

Hi. I am trying to solve a system of differential equations but I am having trouble understanding what I do wrong. I am a novice Maple user.

>with(PDEtools): with(LinearAlgebra): with(DifferentialGeometry): with(VectorCalculus): with(DEtools):

>_Omega_d := 1-_Omega_r

>eq1 := h1(z)*h2(z)+h2(z)*h3(z)+h1(z)*h3(z)-3*(_Omega_d*((1+z)^3)+_Omega_r*((1+z)^4)) = 0

>eq2 := h2(z)^2+h3(z)^2+h2(z)*h3(z)-1/3*(1+z)*(h2(z)*(diff(h2(z), z))+h3(z)*(diff(h3(z), z)))+_Omega_r*((1+z)^4) = 0

>eq3 := h1(z)^2+h3(z)^2+h1(z)*h3(z)-1/3*(1+z)*(h1(z)*(diff(h1(z), z))+h3(z)*(diff(h3(z), z)))+_Omega_r*((1+z)^4) = 0

>eq4 := h2(z)^2+h1(z)^2+h2(z)*h1(z)-1/3*(1+z)*(h2(z)*(diff(h2(z), z))+h1(z)*(diff(h1(z), z)))+_Omega_r*((1+z)^4) = 0

>sys := {eq2, eq3, eq4}

Equations eq1 - eq4 are my equations, but since eq1 does not have a differentiated function, I leave it out of the system, otherwise I get an error that states that I have more functions than variables.

>sys1 := subs(_Omega_r = 10^(-15), sys)

>ics := {h1(0.1e-3) = 1, h2(0.1e-3) = 1, h3(0.1e-3) = 1, (D(h2))(0.1e-3) = 4.500015, (D(h3))(0.1e-3) = 4.500015, (diff(h1(z), z))(0.1e-4) = 4.500015}

ics are my inititial conditions.

>pdsolve(sys1 union ics, {h1(z), h2(z), h3(z)}, type = numeric, generalsolution)

But when I run this, I get an error that says "Error, (in pdsolve/sys) found the independent variables {z} also present in the names of the functions of the system {}".

I have tried numerous different methods to solve this, and I get stuck every time. If anyone has any advice, I would appreciate it. I am trying to get solutions for h1(z), h2(z), h3(z) (to ultimately plot these functions versus z), if possible.