It seems that the analytic result is correct, but the precision is lost when converting it to a number. For example, if we use a higher precision, we get consistent results between numerical and analytical integration:
f[a_, b_] =
Integrate[x^2 Exp[-a x^2 - b x^4], {x, -∞, ∞},
Assumptions -> {a > 0, b > 0}]
g[a_, b_] :=
NIntegrate[x^2 Exp[-a x^2 - b x^4], {x, -∞, ∞}]
ListPlot[{
Table[{a, f[a, 1`50]}, {a, 1/10, 14, 1/10}],
Table[{a, g[a, 1.]}, {a, 1/10, 14, 1/10}]
}, Joined -> True]
You can also set precision use WorkingPrecision in Plot:
Plot[{f[a, 1], g[a, 1]}, {a, 1, 14}, WorkingPrecision -> 50]
