Plotting a small gaussian | Small values and dealing with Machine Precision

I’ve defined the following:

k := 1.38*10^-16 kev := 6.242*10^8 q := 4.8*10^-10 g := 1.66*10^-24 h := 6.63*10^-27 

and

b = ((2^(3/2)) (\[Pi]^2)*1*6*(q^2)*(((1*g*12*g)/(1*g + 12*g))^(   1/2)) )/h  T6 := 20 T := T6*10^6 e0 := ((b k T6 *10^6)/2)^(2/3)  \[CapitalDelta] := 4/\[Sqrt]3 (e0 k T6 *10^6)^(1/2)  \[CapitalDelta]kev = \[CapitalDelta]*kev e0kev = e0*kev bkev = b*kev^(1/2) 

Then, I want to plot these functions:

fexp1[x_] = E^(-bkev *(x*kev)^(-1/2)) fexp2[x_] = E^(-x/(k*T)) fexp3[x_] = fexp1[x]*fexp2[x] 

and check that this Taylor expansion works:

fgauss[x_] =   Exp[(-3 (bkev^2/(4 k T*kev ))^(1/3))]*   Exp[(-((x*kev - e0kev)^2/(\[CapitalDelta]kev/2)^2))] 

which should, e.g., as expected:

Figure 10.1

This plot came from "Stellar Astrophysics notes" of Edward Brown (also it is a known approximation).

I used this line of command to Plot:

Plot[{fexp1[x],fexp2[x],fexp3[x],fgauss[x]}, {x, 0, 50},   PlotStyle -> {{Blue, Dashed}, {Dashed, Green}, {Thick, Red}, {Thick,      Black, Dashed}}, PlotRange -> Automatic, PlotTheme -> "Detailed",   GridLines -> {{{-1, Blue}, 0, 1}, {-1, 0, 1}},   AxesLabel -> {Automatic}, Frame -> True,   FrameLabel -> {Style["EnergĂ­a E", FontSize -> 25, Black],     Style["f(E)", FontSize -> 25, Black]}, ImageSize -> Large,   PlotLegends ->    Placed[LineLegend[{"","","",""}, Background -> Directive[White, Opacity[.9]],      LabelStyle -> {15}, LegendLayout -> {"Column", 1}], {0.35, 0.75}]] 

but it seems that Mathematica doesn’t like huge negative exponentials. I know I can compute this using Python but it’s a surprise to think that Mathematica can’t deal with the problem somehow. Could you help me?