I believe MOL uses the wrong sign when discretizing the first derivate. It should generate the appropriate upwind operator, but gets confused by (xb-x) I reckon.
When I turn the sign around
eq = Dτ(g(τ,x)) ~ 0.5*σ^2*Dxx(g(τ,x)) - λ*(x - xb)*Dx(g(τ,x)) - r*g(τ,x)
the instability is gone (I’ve extended the integration time to 5)
I have also used initial conditions that are commensurable with the boundary conditions.
h(x) = 1.0 + cos(π*(x-x_min)/(x_max-x_min))
Note that these BCs alone do not resolve the instability.