Why is Magnetic Susceptibility not showing the expected transition in 2D&3D Ising Model?

Question

I'm trying to code the Ising Model with the metropolis algorithm to study the ferromagnetic-paramagnetic transitions. The code seems to work ; the equilibration happens. While equilibrating, the graphs for the energy and magnetisation of states look fine. After achieving equilibrium; when calculating the different observables for some range of temperatures; the plot for average Energy and Magnetisation are as expected. The Heat capacity v/s temperature plot is also most of the time as expected. But the Magnetic susceptibility plot always is erroneous; as in, the peak is never on the ~2.2 $k_bt$ value; but rather in a $3-4$$ k_bt$ range.

Till now I've tried the following -

1. Tuned the number of Monte Carlo Steps(MCS) for equilibrating and then calculating. I'm runnng $5\times 10^5$ total MCS and calculating the quantities from the last $10^4$ values only so that I'm sure I'm calculating only over states in equilibrium.
2. Changed the grid size to check whether that's affecting it but it doesn't do anything since I'm running enough MCS to get even the larger grids to equilibrate.
3. Trying to cherry pick and do the calculations on only uncorrelated data by picking every 10th value of the $10^4$ values. I haven't calculated the correlation time so maybe 10 isn't proper, maybe I should try with some more values than 10.

The following is my code-

import numpy as np
import matplotlib.pyplot as plt

N = 50 #sets the grid size 
x = np.random.rand(N,N)


x=np.ones((N,N)) # Generates a array with all +1; i.e , a T=0 state

# Gives us the energy associated with single spin; takes care of periodic boundary conditions
def spin(i,j):
    if i==N-1 and j==N-1:
        return x[i,j]*(x[i,j-1] + x[0,j] + x[i,0] + x[i-1,j]) 
    elif i==N-1:
        return x[i,j]*(x[i,j-1] + x[0,j] + x[i,j+1] + x[i-1,j])
    elif j==N-1:
        return x[i,j]*(x[i,j-1] + x[i+1,j] + x[i,0] + x[i-1,j])
    else:
        return x[i,j]*(x[i,j-1] + x[i+1,j] + x[i,j+1] + x[i-1,j]) 

#Gives us the energy of the state
def Energy(array):
    Enn=0
    for i in range(0,N):
        for j in range(0,N): 
            Enn = Enn - spin(i,j)
    return (1/2)*Enn
    
#Gives us the magnetisation of the state
def mag(array):
    return array.sum()

#This function implements the metropolis algorithm and calculates the needed quantities
def equilibrate(x,T):
    H = Energy(x)
    m = mag(x)
    steps=0
    E_plot=[H]
    magnetisations=[m]
    
    for i in range(5*10**5):
        p = np.random.randint(0,N)
        q = np.random.randint(0,N)
        spin_i=x[p,q]
        steps = steps +1
        dE = 2*spin(p,q)
        if dE <=0:
            x[p,q]=-x[p,q]
            E_new = H + dE
            m_new= m +(-2)*spin_i
        else:
            temp=np.random.rand()
            prob = np.exp(-(dE)/T)
            if temp < prob:
                x[p,q]=-x[p,q]
                E_new = H + dE
                m_new= m +(-2)*spin_i
            else:
                E_new = H
                m_new= m
        H = E_new
        m = m_new
        E_plot.append(E_new)
        magnetisations.append(m_new)

    E_calc = np.array(E_plot)[-10**4:][::100]
    M_calc = np.array(magnetisations)[-10**4:][::100]
    
    E_avg= E_calc.mean()
    C= ((E_calc**2).mean() - E_avg**2)/(T**2 * N**2)
    M_avg= M_calc.mean()
    Chi = ((M_calc**2).mean()-M_avg**2)/T
    
    M_avg = M_avg/ (N**2)
    
    return [E_avg , M_avg, C, Chi]

temp=np.arange(0.1,5,0.1) #Generating a set of temperatures 

E_s=[]
M_s=[]
C_s=[]
Chi_s=[]
for T in temp:
    E_avg,M_avg, C, Chi=equilibrate(x,T)
    E_s.append(E_avg)
    M_s.append(M_avg)
    C_s.append(C)
    Chi_s.append(Chi)
#Running the metropolis algorithm on the state at different temperatures

The code for the 3D case is analogous; just the periodic boundary effects are tweaked a bit and a very similar thing happens there also; though the heat-capacity plots are too somewhat erroneous there.

Here are the plots that I get for 2D -