%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from qutip import *
# Hamiltonian
B = 1.0
H = B * sigmaz()
#initial state
psi0 = (basis(2, 0) + basis(2,1)).unit() # super position : (|0> + |1>)
#psi0 = basis(2, 1)
tlist = np.linspace(0, 10, 100)
#unitary dynamics calculation
result = mesolve(H, psi0, tlist, [], [sigmax(), sigmay(), sigmaz()])
fig, axes = plt.subplots(1,1)
axes.plot(tlist, result.expect[0], label = r'$\left<\sigma_x\right>$')
axes.plot(tlist, result.expect[1], label = r'$\left<\sigma_y\right>$')
axes.plot(tlist, result.expect[2], label = r'$\left<\sigma_z\right>$')
axes.legend(loc=2);