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.