Find best minima for initial conditions for past solution to match present values known

I have in Python an ODE system to solve and it is splitted in 2 parts : a past time solution and a future time solution :

Here is the ODE system :

# defining system
def sys_to_solve(t, funcs):
    return [
        funcs[0] * funcs[1],
        (
            3
            / (3 + 2 * omega_BD)
            * psi_x
            / funcs[0]
            * H0 ** 2
            * Omega_m
            * funcs[2] ** -3
            - (funcs[1]) ** 2
            - 3 * H(funcs[2], funcs[1], funcs[0]) * funcs[1]
        ),
        funcs[2] * H(funcs[2], funcs[1], funcs[0]),
    ]

I know present value (at t_now=0) of the 3 variables.

I am looking for finding initial values ( t=t_begin ) that will match, once solved, the present values.

The known values of present time are x_0, x_1, x_2

Here is the past solution implementation :

sol_past = solve_ivp(
        sys_to_solve,
        [t_begin, t_now],
        [z1, z2 , z3],
        method="RK45",
        t_eval=t_past,
        rtol=1e-12,
        )

To find these initial values, someone told me, to find appropriate z1,z2,z3 above, the following method :

def objective(x):
    sol_past = solve_ivp(
        sys_to_solve,
        [t_begin, t_now],
        [x[0], x[1], x[2]],
        method="RK45",
        t_eval=t_past,
        rtol=1e-12,
    )
    y1 = sol_past.y[0,:]
    y2 = sol_past.y[1,:]
    y3 = sol_past.y[2,:]

    return np.abs(y1[-1] - x_0), np.abs(y2[-1] - x_1), np.abs(y3[-1] - x_2)

z1,z2,z3 = fsolve(objective, [6.45, 1.317, 0.3] , xtol=1e-12)

But the values 6.45, 1.317 and 0.3 have been found by random or step by step methods.

Now, I would like to justify and find these values with a numerical method well established, and not just a pseudo-random fine-tuning.

I have heard this method is called the "shooting method".

I think about minimizers methods like the minimization of a Chi2 : Could anyone help me how to proceed to find similar values 6.45, 1.317 and 0.3 that will match the best the present values x_0,x_1,x_3 ?

EDIT : the solution given by @Lutz Lehmann works well on the 2 first parameters ( Psi = G * phi and F0 = d_Phi_0 / dt / Phi0 ) but I have not much a good math for the scale factor H(t) = da/dt / a given the fact that I have only access to a(t) and not da/dt /a.

Here the bad match on H(t) :

Bad match witht expected H(z) in blue point

However, I get consistent curves for the curve of scale factor and Psi(t) and F0 = d_Phi_0 / Phi0

Well match with a(t

For information, we get funcs[0] = Psi(t), funcs1=F(t) and funcs2=a(t)

So surely, the shooting method could solve this issue on H(t) with currently too high values.

The format is always that the time span is given as [t0,tf] for with initial point y0 for the condition y(t0)=y0 and integration to the time tf.

It makes no difference for the solver if tf > t0 or tf < t0. If given, the t_eval array has to be ordered appropriately. So you only need to change the call to

sol_past = solve_ivp(
        sys_to_solve,
        [t_now, t_begin],
        [z1, z2 , z3],
        method="RK45",
        t_eval=t_past,
        rtol=1e-12,
        )