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?