Summary: ``We present a characteristic initial value approach to calculating the Green function of the Regge-Wheeler and Zerilli equations. We combine well-known numerical methods with newly derived initial data to obtain a scheme which can in principle be generalized to any desired order of convergence. We demonstrate the approach with implementations up to sixth order in the grid spacing. By combining the results of our numerical code with late-time tail expansions and methods of subtracting the direct part of the Green function, we show that the scalar self-force in Schwarzschild spacetime can be computed to better accuracy than previous Green function based approaches. We also demonstrate agreement with frequency domain methods for computing the Green function in the gravitational case. Finally, we apply the Regge-Wheeler and Zerilli Green functions to the computation of the gravitational energy flux.''