I’m using recursion method to do decimation (what I define as s in my program) to a square lattice. I want to show Zeta in function of w.
The problem is, my variable is E* (energy), more specifically w. If I make w vary at the beginning of the program I don’t know how to plot the Zeta. The code is as follow’s:

wmin = -4;
wmax = 4;
s = 5;
n = 2;
\[Eta] = 0.1; t = 1; \[Epsilon] = 0;


U = Table[If[i == j, t, 0], {i, 1, n}, {j, 1, n}];
H = Table[If[i == j - 1, t, 0], {i, 1, n}, {j, 1, n}] + 
   Table[If[i == j + 1, t, 0], {i, 1, n}, {j, 1, n}] +
   Table[If[i == j, \[Epsilon], 0], {i, 1, n}, {j, 1, n}];

This first part, I define some matrices…The loop part to make w vary is:

While[
 wmin <= wmax,

 SubStar[E] = IdentityMatrix[n]*(wmin + \[Eta]*I);
 Subscript[\[Alpha], 1] = U .Inverse[SubStar[E] - H] .U;

 For[j = 2, j < s + 1, j++,
  Subscript[\[Alpha], j] = 
    Subscript[\[Alpha], j - 1].Inverse[SubStar[E] - H - 2 \!\(
\*SubsuperscriptBox[\(\[Sum]\), \(i = 1\), \(j - 1\)]
\*SubscriptBox[\(\[Alpha]\), \(i\)]\)].Subscript[\[Alpha], j - 1];

  ];


 G11 = Inverse[SubStar[E] - H - \!\(
\*SubsuperscriptBox[\(\[Sum]\), \(i = 1\), \(s\)]
\*SubscriptBox[\(\[Alpha]\), \(i\)]\)];

 \[Zeta] = -1/\[Pi] Im[G11];

 wmin = wmin + 1;

 ]

I actually know what I expect from the plot. In other version of the program, I let w vary only at the end of the code (Inside the Plot), and as you might imagine, it’s a veeeery slow way to get a result.

Plot[Diagonal[\[Zeta]], {w, -4, 4}, Frame -> True, 
 Exclusions -> None]

The result from this plot is:

OBS: I have to use a large s to make the graph’s look better. But it’s get way too slowly.

So, if you have any suggestion to plot in a efficient way, I would be thankful!