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