p[m_,v_] := Exp[ -(x-m)^2/(2 v) ] / Sqrt[2 Pi v]
m1 = (1-eps) xb + eps y
m2 = xb
gd = p[m1,1] - (1-eps) p[m2,1]
Limit[ gd/eps , eps -> 0]
gdx = (1 + x y) Exp[-x^2/2]
Plot3D[gdx,{x,-3,3},{y,-3,3}]
gdx1 = ( gdx /. x -> 1 )
gdx2 = ( gdx /. x -> 2 )
Plot[ {gdx1,gdx2},{y,-4,4}]