from qutip import *

# 5-qubit system
N = 5

identity = qeye(2)
X = sigmax()
Y = sigmay()
Z = sigmaz()


def pauli_i(N, i, pauli):
    # pauli: X Y Z
    tmp = [identity for m in range(N)]
    tmp[i]=pauli
    return tensor(tmp)

n=3
paulix_n = pauli_i(N, n, X)

H0 = 0
for n in range(N):
    H0 += pauli_i(N, n, X)
H0 = -0.5*H0

H1 = 0
for n in range(N-1):
    H1 += 0.1 * pauli_i(N, n, Z) * pauli_i(N, n+1, Y)

import numpy as np
def H1_coeff(t, args):
    return 9 * np.exp(-(t / 5.) ** 2)

H = [H0, [H1, H1_coeff]]
t = [0.5, 0.6]

psi0 = tensor([basis(2, 1)] + [basis(2,0) for i in range(N-1)]) # Define initial state
output = mesolve(H, psi0, t)
print(output)