Another way:
eps = 10^-10; end = 10;
sol[n_] := NDSolve[
{g'[r] == a[r]*g[r]/r, a'[r]/r == g[r]^2 - 1, a[eps] == n, g[end] == 1}
,{a[r], g[r]}, {r, eps, end}, Method -> {"Shooting",
"StartingInitialConditions" -> {a[end] == eps, g[end] == (1 - eps)}}];
Plot[Evaluate[g[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end},PlotRange -> All,
PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, PlotLabels -> "g(r)"]
Plot[Evaluate[a[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end},
PlotRange -> All,
PlotLegends -> {"n=1", "n=2", "n=5", "n=10"},
PlotLabels -> "a(r)"]
For larger values e.g.end=20
:
eps = 10^-20; end = 20;
sol[n_] := NDSolve["Code the same",WorkingPrecision -> 25];
Plot["Code the same"] // Quiet(*Error messages form InterpolatingFunction to Quiet *)