I have succesfully plotted the running of the gauge couplings to one loop order with mathematica.

My code is as follows:

`SMbetafunctions = {g1'[t] == 1/(16*Pi^2)*(41/6)*g1[t]^3, g2'[t] == 1/(16*Pi^2)*(-19/6)*g2[t]^3, g3'[t] == 1/(16*Pi^2)*-7*g3[t]^3, g1[0] == 0.35940, g2[0] == 0.64754, g3[0] == 1.1666}; sol = NDSolve[SMbetafunctions, {g1, g2, g3}, {t, 0, 50}]; Plot[Evaluate[{g1[t], g2[t], g3[t]} /. sol], {t, 0, 50}, PlotLegends -> {g1, g2, g3}, PlotRange -> Automatic , PlotTheme -> "Automatic", Frame -> True, FrameLabel -> {"t = ln(μ)", "Coupling Strength"}]] `

The beta functions are defined as

$ \qquad \frac{\partial g_1}{\partial t} = \frac{1}{16\pi}\frac{41}{6}g_1^3$

with $ t = \log(\frac{\mu}{\mu_0})$ .

So am I evaluating the ODE’s over the variable $ t$ . But I want to express the x-axis in the energy $ \mu$ instead.

How can I adjust the x-axis to reflect the energy?

I have tried to rewrite the RGE’s in terms of $ \mu$ instead of $ t$ giving me

$ \frac{\partial g_1}{\partial \mu} = \frac{1}{\mu} \frac{1}{16\pi}\frac{41}{6}g_1^3$

However, when I evaluate this, I get an $ \frac{1}{0}$ error with, e.g., {mu, 1, 10000}.

With python for example I can easily adjust the values of the x-axis by $ \mu = e^{t}\cdot\mu_0$ and which simply translates the values from to energy and reflects it also on the plot.