"Stern-Gerlach 2"

N = 10

a = zero(N)

a[1] = 1 / (3^(2/3) 1.35411793942640046318)
a[2] = -1 / (3^(1/3) 2.67893853470774789827)

for(n,4,N, a[n] = a[n - 3] / (n + 3 - 4) / (n + 2 - 4))

check(infixform(a[4]) == "0.0591713")
check(infixform(a[5]) == "-0.0215683")
check(a[6] == 0)
check(infixform(a[7]) == "0.00197238")
check(infixform(a[8]) == "-0.000513531")
check(a[9] == 0)

Airy = sum(n, 0, N-1, a[n + 1] x^n)

A = (alpha e / (2 hbar))^(1/3)
B = alpha e hbar / (8 m^2)

C = exp(-i e alpha z t / (4 m))
D = exp(-i e B0 t / (2 m))

psi1 = eval(Airy, x, A B t^2 + A z) C D
psi2 = eval(Airy, x, A B t^2 - A z) conj(C) conj(D)

epsilon1 = -hbar^2 / (2 m) d(psi1,z,z) +
           e hbar / (2 m) (B0 + alpha z) psi1 -
           i hbar d(psi1,t)

epsilon2 = -hbar^2 / (2 m) d(psi2,z,z) -
           e hbar / (2 m) (B0 + alpha z) psi2 -
           i hbar d(psi2,t)

-- cancel exponentials

epsilon1 = epsilon1 conj(C) conj(D)
epsilon2 = epsilon2 C D

-- truncate

for(n, N - 2, 2 N,
 epsilon1 = eval(epsilon1, t^n, 0),
 epsilon2 = eval(epsilon2, t^n, 0)
)

for(n, N / 2, N,
 epsilon1 = eval(epsilon1, z^n, 0),
 epsilon2 = eval(epsilon2, z^n, 0)
)

epsilon1
epsilon2