again, I'm not a numerical guy, cuch's the expert, but (and cuch the expert probably knows why this won't work)
I'd use the expression for [$]\frac{\partial L}{\partial h}[$] together with the boundary condition that [$]L=0[$] when [$]h=+\infty[$] and shoot down from [$]h=+\infty[$]
you can probably kludge together an approximation for [$]L[$] for large [$]h[$] and shoot down from that large value of [$]h[$]
[$]\frac{\partial L}{\partial h}\left.\right|_{k,\rho}=-\frac{1}{\sqrt{8\pi}}\exp(-\frac{h^{2}}{2})\;\mathrm{erfc}(\frac{k-\rho h}{\sqrt{2(1-\rho^2)}})[$]
for large [$]h[$], ignore [$]k[$],
[$]\frac{\partial L}{\partial h}\left.\right|_{k,\rho}=-\frac{1}{\sqrt{8\pi}}\exp(-\frac{h^{2}}{2})\;\mathrm{erfc}(-\frac{\rho h}{\sqrt{2(1-\rho^2)}})
=-\frac{1}{\sqrt{8\pi}}\exp(-\frac{h^{2}}{2})\left[2-\mathrm{erfc}(\frac{\rho h}{\sqrt{2(1-\rho^2)}})\right]
[$]
and I've written it like that because there's an expansion of erfc(x) for large x
here
at first cut, replace the erfc by 0 for large [$]h[$] (if you want more terms, use the link above) so for large [$]h[$]
[$]\frac{\partial L}{\partial h}\left.\right|_{k,\rho}=-\frac{1}{\sqrt{2\pi}}\exp(-\frac{h^{2}}{2})[$]
so for large [$]h[$] we can approximate [$]L=\frac{1}{2}\;\mathrm{erfc}(\frac{h}{\sqrt{2}})[$]
if you can follow that, you're smarter than me, but the point was to try and approximate the boundary condition for large [$]h[$], so I'd solve
[$]\frac{\partial L}{\partial h}\left.\right|_{k,\rho}=-\frac{1}{\sqrt{8\pi}}\exp(-\frac{h^{2}}{2})\;\mathrm{erfc}(\frac{k-\rho h}{\sqrt{2(1-\rho^2)}})[$]
shooting down from some large value of [$]h=H[$] with [$]L(H)=\frac{1}{2}\;\mathrm{erfc}(\frac{H}{\sqrt{2}})[$]
It would be even easier to start from [$]L(H)=0[$] but if I've understood, that gives the wrong answers