(*1×1模型，PBCスペクトル*)
Do[t[δx,δy]=RandomReal[{-1,1}],{δx,-1,1},{δy,-1,1}];
h=Sum[E^(-I δx kx-I δy ky)*t[δx,δy],{δx,-1,1},{δy,-1,1}];
Δ=100;
Eval=Table[h,{kx,0,2 Pi,2 Pi/Δ},{ky,0,2 Pi,2 Pi/Δ}]//Flatten;
ListPlot[{Re[#],Im[#]}&/@Eval]

(*2×2模型.PBCスペクトル*)
Do[t[i,j,x,y]=RandomReal[{-1,1}],{i,1,2},{j,1,2},{x,-1,1},{y,-1,1}];
h=Table[Sum[E^(-I x kx-I y ky)*t[i,j,x,y],{x,-1,1},{y,-1,1}],{i,1,2},{j,1,2}];
Δ=100;
Eval=Table[Eigenvalues[h],{kx,0,2 Pi,2 Pi/Δ},{ky,0,2 Pi,2 Pi/Δ}]//Flatten;
ListPlot[{Re[#],Im[#]}&/@Eval]

(*1×1模型．PBCとOBCスペクトル*)
Do[t[δx,δy]=RandomReal[{-1,1}],{δx,-1,1},{δy,-1,1}];
h=Sum[E^(-I δx kx-I δy ky)*t[δx,δy],{δx,-1,1},{δy,-1,1}];
Δ=100;
Eval=Table[h,{kx,0,2 Pi,2 Pi/Δ},{ky,0,2 Pi,2 Pi/Δ}]//Flatten;
ListPlot[{Re[#],Im[#]}&/@Eval]
(*OBC*)
Lx=20;Ly=20;
n[x_,y_]:=Lx (y-1)+x;
hO=Table[0,{i,1,Lx Ly},{j,1,Lx Ly}];
Do[
hO[[n[x,y],n[x+1,y]]]=t[1,0];
,{x,1,Lx-1},{y,1,Ly}];
Do[
hO[[n[x,y],n[x+1,y+1]]]=t[1,1];
,{x,1,Lx-1},{y,1,Ly-1}];
Do[
hO[[n[x,y],n[x+1,y-1]]]=t[1,-1];
,{x,1,Lx-1},{y,2,Ly}];
Do[
hO[[n[x,y],n[x,y]]]=t[0,0];
,{x,1,Lx},{y,1,Ly}];
Do[
hO[[n[x,y],n[x,y+1]]]=t[0,1];
,{x,1,Lx},{y,1,Ly-1}];
Do[
hO[[n[x,y],n[x,y-1]]]=t[0,-1];
,{x,1,Lx},{y,2,Ly}];
Do[
hO[[n[x,y],n[x-1,y]]]=t[-1,0];
,{x,2,Lx},{y,1,Ly}];
Do[
hO[[n[x,y],n[x-1,y+1]]]=t[-1,1];
,{x,2,Lx},{y,1,Ly-1}];
Do[
hO[[n[x,y],n[x-1,y-1]]]=t[-1,-1];
,{x,2,Lx},{y,2,Ly}];

Eval=Eigenvalues[hO];
ListPlot[{Re[#],Im[#]}&/@Eval]
(*全固有関数の密度プロット*)
Efunc=Eigensystem[hO][[2]];
Density=Table[Total[Norm/@Efunc[[All,n[x,y]]]]/(Lx*Ly),{x,1,Lx},{y,1,Lx}];
ListPlot3D[Flatten[Table[{x,y,Density[[x,y]]},{x,1,Lx},{y,1,Ly}],1]]