diff --git a/presentation/spin_chain/hamiltonian.py b/presentation/spin_chain/hamiltonian.py
index c33a6c6..8274bee 100644
--- a/presentation/spin_chain/hamiltonian.py
+++ b/presentation/spin_chain/hamiltonian.py
@@ -11,6 +11,7 @@ def Mi(nqbits, i, M):
             result = np.kron(result, I)
         else:
             result = np.kron(result, M)
+    return result
 
 
 
@@ -23,5 +24,5 @@ def H_field(nqbits, g):
     return sum(field_terms)
 
 def H(nqbits, g):
-    return (-H_interaction + H_field).real
+    return (-H_interaction(nqbits) + H_field(nqbits, g)).real
 
diff --git a/presentation/spin_chain/time_evolution.py b/presentation/spin_chain/time_evolution.py
index 9c4f639..712d65f 100644
--- a/presentation/spin_chain/time_evolution.py
+++ b/presentation/spin_chain/time_evolution.py
@@ -5,6 +5,8 @@ from pyqcs import State, sample
 from transfer_matrix import T_time_slice
 from hamiltonian import H
 
+from scipy.linalg import expm
+
 nqbits = 4
 g = 0.5
 N = 400
@@ -17,6 +19,7 @@ measure = 0b10
 
 
 results_qc = []
+results_np = []
 print()
 for t in np.arange(0, t_stop, delta_t):
     # QC simulation
@@ -30,14 +33,27 @@ for t in np.arange(0, t_stop, delta_t):
     results_qc.append(result[0] / n_sample)
 
     # Simulation using matrices
-    #np_
+    np_zero_state = np.zeros(2**nqbits)
+    np_zero_state[0] = 1
+
+    T = expm(-1j * t * H(nqbits, g))
+    #for Tv in T:
+    #    print(np.sum(np.abs(Tv)))
+    #    assert np.isclose(np.sum(np.abs(Tv)), 1)
+    np_state = T.dot(np_zero_state)
+    amplitude = np.sum(np.abs(np_state[[False if (i & measure) else True for i in range(2**nqbits)]]))
+    results_np.append(amplitude)
+
     print(f"simulating... {int(t/t_stop*100)} %  ", end="\r")
 print()
 print("done.")
 
-plt.plot(np.arange(0, t_stop, delta_t), results_qc)
+h0, = plt.plot(np.arange(0, t_stop, delta_t), results_qc, label=f"Quantum computing {n_sample} samples")
+h1, = plt.plot(np.arange(0, t_stop, delta_t), results_np, label="Classical simulation using explicit transfer matrix")
 plt.xlabel("t")
 plt.ylabel(r"$|0\rangle$ probability amplitude for second spin")
 plt.title(f"{nqbits} site spin chain with g={g} coupling to external field")
+plt.legend(handles=[h0, h1])
+
 plt.show()