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:

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?