I want to calculate Coulomb integrals with a double continuum in HF. I know that this method is actually only for bound states. But is there any possibility to approximate some wave functions in the continuum and adding them to my basic set??

The long answer is really long, because there are various ways to do this, from the poor man’s approach to very fancy.
Very fancy: search complex scaling or complex absorbing potential.

More straightforward, the Hazi-Taylor stabilization method. Effectively, any specific calculation models a discretized contimuum, that is a pseudo continuum modeled by “particle in the box” states, where the basis set defines an implicit box. What you do then is to systematically vary the size of this basis set box, and to average over it, and that will give you the right answer for the continuum. Check: @article{mandelshtam93,
author = “V. A. Mandelshtam and T. R. Ravuri and H. S. Taylor”,
journal = prl,
volume = 70,
pages = 1932,
year = 1993,
title = “Calculation of the density of resonance states using the stabilization method”}
and @article{mandelshtam94,
author = “V. A. Mandelshtam and T. R. Ravuri and H. S. Taylor”,
journal = jcp,
volume = 101,
pages = 8792,
year = 1994,
title = “The stabilisation theory of scattering”}

The poor man’s implementation of this to vary the exponents of 2 or 3 extra diffuse basis functions. Personally, I do not like this approach much, because it changes the basis sets. I’s rather keep my basis set the same and add an explicit confining potential to the Hamiltonian, which is not too much effort in electron-scattering at the HF level, but, of course, less trivial at any many-body level.