text: import numpy as np
  import matplotlib.pyplot as plt
  import time
  
  from numba import jit
  
  
  @jit(nopython=True)
  def energy(system, i, j, L):
      """Energy function of spins connected to site (i, j)."""
      return -1. * system[i, j] * (system[np.mod(i - 1, L), j] + system[np.mod(i + 1, L), j] +
                                   system[i, np.mod(j - 1, L)] + system[i, np.mod(j + 1, L)])
  
  
  @jit
  def prepare_system(L):
      """Initialize the system."""
      system = 2 * (0.5 - np.random.randint(0, 2, size=[L, L]))
      return system
  
  @jit(nopython=True)
  def measure_energy(system):
      L = system.shape[0]
      E = 0
      for i in range(L):
          for j in range(L):
              E += energy(system, i, j, L) / 2.
      return E
  
  
  @jit(nopython=True)
  def metropolis_loop(system, T, N_sweeps, N_eq, N_flips):
      """ Main loop doing the Metropolis algorithm."""
      E = measure_energy(system)
      L = system.shape[0]
      E_list = []
      for step in range(N_sweeps + N_eq):
          i = np.random.randint(0, L)
          j = np.random.randint(0, L)
  
          dE = -2. * energy(system, i, j, L)
          if dE <= 0.:
              system[i, j] *= -1
              E += dE
          elif np.exp(-1. / T * dE) > np.random.rand():
              system[i, j] *= -1
              E += dE
  
          if step >= N_eq and np.mod(step, N_flips) == 0:
              # measurement
              E_list.append(E)
      return np.array(E_list)
  
  
  
  
  if __name__ == "__main__":
      """ Scan through some temperatures """
      # Set parameters here
      L = 4  # Linear system size
      N_sweeps = 5000  # Number of steps for the measurements
      N_eq = 1000  # Number of equilibration steps before the measurements start
      N_flips = 10  # Number of steps between measurements
      N_bins = 10  # Number of bins use for the error analysis
  
      T_range = np.arange(1.5, 3.1, 0.1)
  
      C_list = []
      system = prepare_system(L)
      for T in T_range:
          C_list_bin = []
          for k in range(N_bins):
              Es = metropolis_loop(system, T, N_sweeps, N_eq, N_flips)
  
              mean_E = np.mean(Es)
              mean_E2 = np.mean(Es**2)
  
              C_list_bin.append(1. / T**2. / L**2. * (mean_E2 - mean_E**2))
          C_list.append([np.mean(C_list_bin), np.std(C_list_bin) / np.sqrt(N_bins)])
  
          print(T, mean_E, C_list[-1])
  
      # Plot the results
      C_list = np.array(C_list)
      plt.errorbar(T_range, C_list[:, 0], C_list[:, 1])
      Tc = 2. / np.log(1. + np.sqrt(2))
      print(Tc)
      plt.axvline(Tc, color='r', linestyle='--')
      plt.xlabel('$T$')
      plt.ylabel('$c$')
      plt.show()
  