Strange phenomenon occurring in analytic integration result involving Bessel functions

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 the precision in Plot using WorkingPrecision:

Plot[{f[a, 1], g[a, 1]}, {a, 1, 14}, WorkingPrecision -> 50]

The culprit, as suspected by xslittlegrass, is indeed numerical instability; in particular, this is because of the perverse combination of modified Bessel functions exhibited in the result returned by Mathematica.

Using a recurrence identity satisfied by the modified Bessel function of the first kind, we can simplify the expression returned, like so:

Integrate[x^2 Exp[-a x^2 - b x^4], {x, -∞, ∞}, Assumptions -> {a > 0, b > 0}] /.
          BesselI[5/4, a^2/(8 b)] ->
          BesselI[-3/4, a^2/(8 b)] - 4 b/a^2 BesselI[1/4, a^2/(8 b)] // FullSimplify
   (a^2 E^(a^2/(8 b)) (- BesselK[1/4, a^2/(8 b)] + BesselK[3/4, a^2/(8 b)]))/(8 Sqrt[a b^3])

where the result is now in terms of the function of the second kind.

(As to why Mathematica is unable to return this result automatically, I've no idea.)

Thus,

pcefunction[a_, b_] :=
   (a^2 E^(a^2/(8 b)) (BesselK[3/4, a^2/(8 b)] - BesselK[1/4, a^2/(8 b)]))/(8 Sqrt[a b^3])

Plot[pcefunction[a, 1], {a, 0, 15}]