This can indeed be done with the Polya Enumeration Theorem. There is a very simple yet powerful observation that can be made here: if we have the cycle index of the edge permutation group induced by the action of the symmetric group on the vertices, the cycle index of the edge permutation group including self-loops can be obtained by multiplying the factorizations into disjoint cycles of the edge permutation and the original vertex permutation and sum these contributions for all vertex permutations. (Of course we work directly with the cycle index of the symmetric group which is much more efficient than iterating over all permutations.) But the cycle index of the edge permutation group is not difficult to compute, the code can be found at this MSE Link.
This gives the following cycle indices for the enumeration of graphs that may have self-loops: for $n=2$ we have $1/2\,{a_{{1}}}^{3}+1/2\,a_{{1}}a_{{2}},$ for $n=3$ we get $1/6\,{a_{{1}}}^{6}+1/2\,{a_{{1}}}^{2}{a_{{2}}}^{2}+1/3\,{a_{{3}}}^{2},$ and finally for $n=4$ we obtain $1/24\,{a_{{1}}}^{10}+1/4\,{a_{{1}}}^{4}{a_{{2}}}^{3}+1/3\,{a_{{3}}}^{3}a_{{1}}+1 /8\,{a_{{1}}}^{2}{a_{{2}}}^{4}+1/4\,a_{{2}}{a_{{4}}}^{2}.$
As an example we substitute $1+z$ into the cycle index for $n=3$ to get ${z}^{6}+2\,{z}^{5}+4\,{z}^{4}+6\,{z}^{3}+4\,{z}^{2}+2\,z+1.$ There are indeed two graphs on three vertices with $1$ edge including loops, namely a single edge or a single loop. Similarly there are four graphs with $2$ edges, namely two loops, two edges, an edge with a loop on one end and an edge that is opposed to a vertex with a loop on it. Furthermore the six graphs with three edges are: a triangle with no loops, a fork with a loop on the handle, a fork with a loop on one of its ends, a line segment with two loops at either end, a line segment with a loop on one end and a loop not incident on the line segment, and finally, three loops. It goes on like this.
These cycle indices produce the following total counts of non-isomorphic undirected graphs with possible self-loops present:
1 6 20 90 544 5096 79264 2208612 113743760 10926227136 1956363435360 652335084592096 405402273420996800 470568642161119963904 1023063423471189431054720 4178849203082023236058229792 32168008290073542372004082199424 468053896898117580623237189908861696 12908865186281154904787051023018388368640 676599895346467962508983189158040730734206464 67557268911383303274368343026941404659790140175872 12878644470123443279570180725329554086442149175832124416 4696891990987444795146988875290693218960182452250871238964224 3283275444870880220030739747435965751837707129958501790119044590592 4406671511641658245922265014648856899986657018190959680004287879813376000 11374011362523188765310166602160674880959112891073505609667787021125321629659136
This sequence of values points us to OEIS A000666 where a large amount of additional material can be found.
For those who are curious to see the details, here is the Maple code to compute these values:
with(numtheory); with(group): with(combinat): pet_cycleind_vrec := proc(n) local p, s; option remember; if n=0 then return 1; fi; expand(1/n*add(a[l]*pet_cycleind_vrec(n-l), l=1..n)); end; pet_flatten_term := proc(varp) local terml, d, cf, v; terml := []; cf := varp; for v in indets(varp) do d := degree(varp, v); terml := [op(terml), seq(v, k=1..d)]; cf := cf/v^d; od; [cf, terml]; end; pet_cycleind_edg := proc(n) option remember; local s, t, res, cycs, l1, l2, flat, u, v; if n=0 then return 1; fi; s := 0: for t in pet_cycleind_vrec(n) do flat := pet_flatten_term(t); cycs := flat[2]; res := 1; for u to nops(cycs) do for v from u+1 to nops(cycs) do l1 := op(1, cycs[u]); l2 := op(1, cycs[v]); res := res * a[lcm(l1, l2)]^(l1*l2/lcm(l1, l2)); od; od; for u to nops(cycs) do l1 := op(1, cycs[u]); if type(l1, odd) then # a[l1]^(1/2*l1*(l1-1)/l1); res := res*a[l1]^(1/2*(l1-1)); else # a[l1/2]^(l1/2/(l1/2))*a[l1]^(1/2*l1*(l1-2)/l1) res := res*a[l1/2]*a[l1]^(1/2*(l1-2)); fi; od; s := s + res*t; od; s; end; pet_varinto_cind := proc(poly, ind) local subs1, subs2, polyvars, indvars, v, pot, res; res := ind; polyvars := indets(poly); indvars := indets(ind); for v in indvars do pot := op(1, v); subs1 := [seq(polyvars[k]=polyvars[k]^pot, k=1..nops(polyvars))]; subs2 := [v=subs(subs1, poly)]; res := subs(subs2, res); od; res; end; v := proc(n) option remember; local p, k, gf; p := q1+q2; gf := expand(pet_varinto_cind(p, pet_cycleind_edg(n))); subs({q1=1, q2=1}, gf); end;