quantum mechanics - Numerical approximation of the wavefunction in a delta-potential

I am trying to approximate the wavefunction of a particle in a delta potential $U(x) = -U_0 \delta(x)$ with $V_0 \gt 0$. I am using the following formula to calculate the wavefunction:

$\psi(x+\Delta x) = 2 \psi(x) - \psi(x-\Delta x) + (\Delta x)^2 (V(x) - \epsilon)\psi(x)$

This can be derived from the one dimensional, time-independent, one-dimensional schrödinger-equation $\psi''(x) = (U(x)-\epsilon)\psi(x).$

The following C++-code sippet shows the calculation:

    double phi_x_prev = 0;
    double phi_x_curr = 1e-100;
    double phi_x_next = 0;

    m_value.push_back(phi_x_curr);

    double stepSq = step * step;

    for (double x = -xmax + step; x <= xmax; x += step)

    {
        phi_x_next =
            (2 * phi_x_curr - phi_x_prev)                           // O(1)
            + stepSq * (system->Potential(x) - eps)*phi_x_curr;     // O(step^2)

        m_value.push_back(phi_x_next);

        phi_x_prev = phi_x_curr;
        phi_x_curr = phi_x_next;
    }

The solution looks good for $x<0$, but for $x>0$ the wavefunction starts growing exponentially, which is not compatible with the boundary condition $\psi(x) \rightarrow \infty$ for $x \rightarrow \pm \infty$.

What did I do wrong? Thanks in advance!