How to calculate initial conditions to integrate a null geodesic Suppose, this is the line element of a FLRW metric,
 $$ ds^2 = -[1 + 2ψ(t,x_i)]dt^2 + a^2(t) [1 - 2ϕ(t,x_i)]dx_i^2 $$
and the geodesic equation is,
$$ \frac{d^2x^α}{dλ^2} = - Γ_{βγ}^α \frac{dx^β}{dλ} \frac{dx^γ}{dλ} $$ 
where λ is the affine parameter. Again, for null geodesic we can write
$$ ds^2 = 0 $$ 
Now, if I want to solve this null geodesic equation I know that I need to convert the four 2nd order differential equation into eight first order equation first and then easily anyone can solve this by using any mathematical tools. Actually, my concern is about the initial conditions, how can I set initial conditions to solve this null geodesic equation, i.e. what will be the initial conditions here for x0, y0, z0 and  for the four velocities,  \$$ \frac{dt0}{dλ}, \frac{dx0}{dλ}, \frac{dy0}{dλ}, \frac{dz0}{dλ}$ respectively, suppose, the initial condition for time is today that means I want to integrate the geodesic equation from today i.e, the redshift is zero. 
 A: There are two parts to your question, let's call them the math and the physics. In what follows a dot will mean a derivative with respect to $\lambda$.
Math
You have, in principle, eight numbers to specify: $(t_0, x_0, y_0, z_0, \dot{t}_0, \dot{x}_0, \dot{y}_0, \dot{z}_0)$. The initial positions and time are just four numbers you have to choose, but the velocities are a bit more subtle. This is because the four-velocity at all times must satisfy $g_{\mu\nu} \dot{x}^\mu \dot{x}^\nu = 0$. This lets you solve for one component in terms of the other three: we usually specify the spatial velocities and solve for $\dot{t}_0$, since this goes with our intuitive idea of saying the direction in which the particle is moving.
With null geodesics, however, there is one more step. Since we can reparametrize $\lambda \to a \lambda + b$ while keeping the null condition $g_{\mu\nu} \dot{x}^\mu \dot{x}^\nu = 0$ invariant, multiplying the whole vector $\dot{x}^\mu$ by a constant gives us a physically equivalent vector. This means that we really have two degrees of freedom: one of the $\dot{x}^i_0$ can be arbitrarily fixed, and only the other two have physical meaning. One simple way to deal with this is to go to spherical coordinates. There are a few ways to do it, but the easiest is just to pretend that $\dot{x}^i_0$ are the components of a regular 3D Euclidean vector and use the usual conversion. That way, $\theta$ and $\varphi$ are the physically meaningful quantities, and $r$ is irrelevant (set it to $1$ or whatever).
To sum up, a null geodesic really has only six degrees of freedom: the four components of space-time position, and two numbers specifying the direction of the three-velocity. Symmetries can make some of these irrelevant: for a example, in a stationary space-time all values of $t_0$ give the same curve.
Physics
You might still be wondering, but what initial conditions do I set? The problem is that there's no such thing as the initial conditions. This has nothing to do with GR. Suppose, for example, that you have an old wooden floor. Moisture has gotten everywhere, so that the floor is all bumpy, higher in some places and lower in others. Now you ask: "if I roll a ball on this floor, what will its path be?". Well, we can't answer that. We're missing information: where you throw the ball from, and in which direction. The ball can take all kinds of possible paths.
The same is true in your spacetime. It's neither homogeneous nor isotropic, so the path of a photon depends on where (and when) you shot it from, in addition to its direction. You might be getting confused because the examples we first learn have a lot of symmetry. If we take for example the Schwarzschild solution, which is stationary and spherically symmetric, we can eliminate a lot of degrees of freedom, and the result is that a photon really only has two: the initial radius and the angular momentum (or impact parameter). Changing any of the other initial conditions just gives a rotated and/or time translated version of the trajectory. This is not the case here; without symmetries, all six numbers are necessary.
The only answer I can give is that you need to define better why you want to integrate the equations. Are you trying to get a feel for possible photon paths? Then your best bet is to pick a few different values and see what happens.
