There are two generic sources of errors.
The first is the choice of guess energy for your solution.
Quantization of energy is the result of applying boundary conditions, and it is by forcing the boundary conditions that one recovers discrete energy levels. If you do this in the reverse direction, i.e. you start with some energy and solve the Schrodinger equation, there is no guarantee that the boundary condition (here at $\infty$) will be satisfied unless you precisely pick one of the solution energies: for pretty much any other value of energy the solution will eventually diverge. The tolerance on $\epsilon$ is quite small: you can be off by less than $1\%$ and the solution will eventually diverge.
The second is numerical accuracy of the scheme. Even if you start with the exact and correct value of energy, the integration scheme will pick up some small round off errors as you move away from your initial conditions. This produces the same effect as picking a slightly wrong value of energy: the solution will behave for a while but eventually, for some large enough values of the argument, will start to exponentially diverge. There is no solution to this beyond using increasingly accurate integration scheme, either by increasing the order of the scheme, or by making the integration mesh finer.
As examples of these two types of error, consider the numerical solution to the Schrodinger equation for a harmonic potential. On the first figure, I solved the equation using the exact eigenvalue for $n=2$: the solution still eventually diverges; this illustrates the second type of error. On the bottom I took $n=2.001$, i.e. an error in guess eigenvalue of 1 part in 2000. Note how the solution starts to diverge for smaller $\vert x\vert$. This shows how a generic integration scheme is very sensitive to the guess eigenvalue; this illustrates the first source of errors. Both integrations were done using the same
NDSolve function in Mathematica so there can be no issue of implementation.
In your specific potential the control of numerical errors is critical because the derivative is discontinuous at the $\delta$. I suspect that some small numerical error occur when going across the $\delta$ and this messes up your solution when you continue to integrate on one side.