How to simulate in python the superconducting density of states of a BCS superconductor with YSR (yu-shiba-rusinov) with 1D chain tight binding model

Question

I try to simulate the following problem:

Compute the density of states numerically of a BCS superconductor with magnetic impurities using the tight binding model of 1D chain with periodically border conditions.

In the image there is a picture of the Hamiltonian. ​ Hamiltonian and matrix to diagonalize

To solve it I want to write the matrix, diagonalize it and use the solutions of eigenvalues and eigenvectors to compute the density of states

Density of states formula

This is the code I have:

import numpy as np
import matplotlib.pyplot as plt

#parameters
N = 2000  # matrix size
delta = 1  # gap
t = 5.0  #  hopping
U = 20.0  
JS = 17  
eta = 0.01 # broadening

# initialize the matrix with zeros
H = np.zeros((N, N), dtype=complex)

# Fill the matrix with the gap and hoping values
for i in range(0, N - 4, 4):
    H[i, i + 3] = -delta
    H[i + 1, i + 2] = delta
    H[i + 2, i + 1] = np.conj(delta)
    H[i + 3, i] = -np.conj(delta)
    H[i, i + 4] = -t
    H[i + 1, i + 5] = -t
    H[i + 2, i + 6] = t
    H[i + 3, i + 7] = t
    H[i + 4, i] = -t
    H[i + 5, i + 1] = -t
    H[i + 6, i + 2] = t
    H[i + 7, i + 3] = t


H[N-4, 0] = -t
H[N-3, 1] = -t
H[N-2, 2] = t
H[N-1, 3] = t
H[0, N-4] = -t
H[1, N-3] = -t
H[2, N-2] = t
H[3, N-1] = t

H[0, 0] = U + JS
H[1, 1] = U - JS
H[2, 2] = -(U + JS)
H[3, 3] = -(U - JS)

# diagonalize the matrix
eigenvalues, eigenvectors = np.linalg.eigh(H)

# Compute the LDOS in one site
def compute_ldos(E_vals, E, psi, eta, site):
    """ Calcula la LDOS en un sitio específico. """
    a_i = psi[4 * site, :]
    b_i = psi[4 * site + 1, :]
    c_i = psi[4 * site + 2, :]
    d_i = psi[4 * site + 3, :]
    LDOS = (np.abs(a_i)**2 + np.abs(b_i)**2) * (eta / np.pi) / ((E - E_vals)**2 + eta**2) + \
           (np.abs(c_i)**2 + np.abs(d_i)**2) * (eta / np.pi) / ((E + E_vals)**2 + eta**2)
    return np.sum(LDOS)


# 
E_range = np.linspace(-3, 3 , 2000)

# Calculate the density of states in site 250
site_index =250
LDOS_site = np.array([compute_ldos(eigenvalues, E, eigenvectors, eta, site_index) for E in E_range])


plt.figure(figsize=(8, 6))
plt.plot(E_range, LDOS_site, label=f'Densidad de estados en sitio {site_index}')
plt.xlabel('E')
plt.ylabel('LDOS')

plt.legend()
plt.grid()
plt.xlim(-3,3)

plt.show()

This is what I get when I try to compute the density of states in site 250:

Clearly my code doesn't work as expected, I would expect it to obtain something close to BCS LDOS. I don't understand where the oscillations are coming from and why far away from the gap it doesn't normalize to one.

My question is where are the mistakes in my code? Also, I saw that with N=8000, it takes too long to compute, so I also wonder if there is a more efficient way to it.