Your lattice consists of 5 x 5 x 5 = 125 spins, so your number of Montecarlo steps to reach equilibrium should be >> 125, because you randomly picking a site and flipping it, so random numbers should uniformly generated so that it will cover whole lattice. For much finer measurement of thermodynamic quantities, you should take more number of points between the temperature range and more number of Montecarlo steps >> 25.

So I modified your code a bit, instead of picking a site randomly and flipping, I swept through each site starting row and column wise.

def mcmove(spin_states, beta):
    '''Monte Carlo move using Metropolis algorithm '''
    for x in range(len(spin_states)):
        for y in range(len(spin_states)):
            for z in range(len(spin_states)):
                s = spin_states[x,y,z]
                cost = 2*s*find_neighbours(spin_states,N,x,y,z)
                if cost < 0:
                    s *= -1
                elif rand() < np.exp(-cost*beta):
                    s *= -1
                spin_states[x, y,z] = s
    return spin_states

I ran the simulations on a 16 x 16 x 16 lattice, temperature ranging from 0.1 - 8.0 with 500 points in-between, with 5000 steps to reach equilibrium and 5000 steps to take measurement. These are my results

Your lattice consists of 5 x 5 x 5 = 125 spins, so your number of Montecarlo steps to reach equilibrium should be >> 125, because you randomly picking a site and flipping it, so random numbers should uniformly generated so that it will cover whole lattice. For much finer measurement of thermodynamic quantities, you should take more number of points between the temperature range and more number of Montecarlo steps >> 25.

So I modified your code a bit, instead of picking a site randomly and flipping, I swept through each site starting row and column wise.

def mcmove(spin_states, beta):
    '''Monte Carlo move using Metropolis algorithm '''
    for x in range(len(spin_states)):
        for y in range(len(spin_states)):
            for z in range(len(spin_states)):
                s = spin_states[x,y,z]
                cost = 2*s*find_neighbours(spin_states,N,x,y,z)
                if cost < 0:
                    s *= -1
                elif rand() < np.exp(-cost*beta):
                    s *= -1
                spin_states[x, y,z] = s
    return spin_states

I ran the simulations on a 16 x 16 x 16 lattice, temperature ranging from 0.1 - 8.0 with 500 points in-between, with 5000 steps to reach equilibrium and 5000 steps to take measurement (you can get similar results with less number of Montecarlo steps, say 100 + 100, since it will sweep through whole lattice). These are my results