diff --git a/presentation/spin_chain/hamiltonian.py b/presentation/spin_chain/hamiltonian.py
new file mode 100644
index 0000000..c33a6c6
--- /dev/null
+++ b/presentation/spin_chain/hamiltonian.py
@@ -0,0 +1,27 @@
+import numpy as np
+
+I = np.matrix([[1, 0], [0, 1]])
+Z = np.matrix([[1, 0], [0, -1]])
+X = np.matrix([[0, 1], [1, 0]])
+
+def Mi(nqbits, i, M):
+    result = 1
+    for j in range(nqbits):
+        if(j != i):
+            result = np.kron(result, I)
+        else:
+            result = np.kron(result, M)
+
+
+
+def H_interaction(nqbits):
+    interaction_terms = [Mi(nqbits, i, Z) @ Mi(nqbits, i+1, Z) for i in range(nqbits)]
+    return sum(interaction_terms)
+
+def H_field(nqbits, g):
+    field_terms = [g*Mi(nqbits, i, X) for i in range(nqbits)]
+    return sum(field_terms)
+
+def H(nqbits, g):
+    return (-H_interaction + H_field).real
+
diff --git a/presentation/spin_chain/time_evolution.py b/presentation/spin_chain/time_evolution.py
new file mode 100644
index 0000000..9c4f639
--- /dev/null
+++ b/presentation/spin_chain/time_evolution.py
@@ -0,0 +1,43 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from pyqcs import State, sample
+from transfer_matrix import T_time_slice
+from hamiltonian import H
+
+nqbits = 4
+g = 0.5
+N = 400
+t_stop = 9
+delta_t = 0.05
+qbits = list(range(nqbits))
+
+n_sample = 200
+measure = 0b10
+
+
+results_qc = []
+print()
+for t in np.arange(0, t_stop, delta_t):
+    # QC simulation
+    state = State.new_zero_state(nqbits)
+
+    for _ in range(N):
+        state = T_time_slice(qbits, t, g, N) * state
+
+    result = sample(state, measure, n_sample)
+
+    results_qc.append(result[0] / n_sample)
+
+    # Simulation using matrices
+    #np_
+    print(f"simulating... {int(t/t_stop*100)} %  ", end="\r")
+print()
+print("done.")
+
+plt.plot(np.arange(0, t_stop, delta_t), results_qc)
+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.show()
+
diff --git a/presentation/spin_chain/transfer_matrix.py b/presentation/spin_chain/transfer_matrix.py
new file mode 100644
index 0000000..0b060e3
--- /dev/null
+++ b/presentation/spin_chain/transfer_matrix.py
@@ -0,0 +1,23 @@
+from pyqcs import X, Z, H, R, CX, State, list_to_circuit
+
+def T_interaction(a, b, t):
+    theta = -t/2
+
+    return (CX(a, b) | R(a, -theta)
+            | X(a) | R(a, theta) | X(a) | CX(a, b))
+
+def T_field(a, t, g):
+    theta = g*t/2
+
+    return (H(a) | R(a, -2*theta) | H(a)
+            | R(a, theta) | X(a) | R(a, theta) | X(a))
+
+def T_time_slice(qbits, t, g, N):
+    interactions_half = list_to_circuit(
+                [T_interaction(i, i+1, t/(2*N))
+                    for i,_ in enumerate(qbits[:-1])]
+            )
+
+    field = list_to_circuit([T_field(i, t/N, g) for i,_ in enumerate(qbits)])
+
+    return (interactions_half | field | interactions_half)