For non-autonomous difference equations of the form x n+1 = f (x n , λ n ), n ∈ we consider homoclinic trajectories. These are pairs of trajectories that converge in both time directions towards each other. Assuming hyperbolicity, we derive a numerical method to compute homoclinic trajectories in two steps. In the first step one trajectory is approximated by the solution of a boundary value problem and precise error estimates are given. In particular, influences of parameters λ n with |n| large are discussed in detail. A second trajectory that is homoclinic to the first one is computed in a subsequent step as follows. We transform the original system into a topologically equivalent form having zero as an n-independent fixed point. Applying the boundary value ansatz to the transformed system, we obtain a non-autonomous homoclinic orbit, converging towards the origin, cf. Hüls (2006). Transforming back to the original coordinates leads to the desired homoclinic trajectories. The numerical method and the validity of the error estimates are illustrated by examples.