This script can be pasted into the Eigenmath app.

 -- Feynman and Hibbs problem 4-2 e = quote(e) -- clear e (was natural number e = exp(1)) A = (A1(),A2(),A3()) x = (x1,x2,x3) eta = (eta1,eta2,eta3) -- Directional derivative (unnormalized) D(f) = dot(eta,d(f,x)) -- Gaussian integral G(f) = do( f = eval(f, eta1^2, i hbar epsilon / m), f = eval(f, eta2^2, i hbar epsilon / m), f = eval(f, eta3^2, i hbar epsilon / m), f = eval(f, eta1, 0), f = eval(f, eta2, 0), f = eval(f, eta3, 0), f (2 pi i hbar epsilon / m)^(3/2)) -- return value -- Factors in integrand u1 = 1 u2 = -i e / (hbar c) dot(eta,A) u3 = -i e / (2 hbar c) dot(eta,D(A)) u4 = 1/2 (-i e / (hbar c) dot(eta,A))^2 u = (u1,u2,u3,u4) v1 = psi() v2 = D(psi()) v3 = 1/2 D(D(psi())) v = (v1,v2,v3) okay(n) = test(n == 1, "ok", "?") "Verify integral 1" f = u[1] v[1] B = G(f) C = psi() (2 pi i hbar epsilon / m)^(3/2) okay(B == C) "Verify integral 2" f = u[1] v[2] B = G(f) okay(B == 0) "Verify integral 3" f = u[1] v[3] B = G(f) C = i hbar epsilon / (2 m) * (2 pi i hbar epsilon / m)^(3/2) * (d(psi(),x1,x1) + d(psi(),x2,x2) + d(psi(),x3,x3)) okay(B == C) "Verify integral 4" f = u[2] v[1] B = G(f) okay(B == 0) "Verify integral 5" f = u[2] v[2] B = G(f) C = -i hbar epsilon / m * (2 pi i hbar epsilon / m)^(3/2) * i e / (hbar c) * dot(A,d(psi(),x)) okay(B == C) "Verify integral 6" f = u[2] v[3] B = G(f) okay(B == 0) "Verify integral 7" f = u[3] v[1] B = G(f) C = -i hbar epsilon / m * (2 pi i hbar epsilon / m)^(3/2) * i e / (2 hbar c) * (d(A1(),x1) + d(A2(),x2) + d(A3(),x3)) * psi() okay(B == C) "Verify integral 8" f = u[3] v[2] B = G(f) okay(B == 0) "Verify integral 9" f = u[3] v[3] B = G(f) B = eval(B,epsilon^(7/2),0) -- discard epsilon^2 okay(B == 0) "Verify integral 10" f = u[4] v[1] B = G(f) C = i hbar epsilon / m * (2 pi i hbar epsilon / m)^(3/2) * 1/2 (i e / (hbar c))^2 * dot(A,A) * psi() okay(B == C) "Verify integral 11" f = u[4] v[2] B = G(f) okay(B == 0) "Verify integral 12" f = u[4] v[3] B = G(f) B = eval(B,epsilon^(7/2),0) -- discard epsilon^2 okay(B == 0) "Verify with Feynman and Hibbs equation (4.18)" -- sum over Gaussian integrals for all uv B = sum(j,1,4, sum(k,1,3, G( u[j] v[k] ) )) -- epsilon^2 epsilon^(3/2) == epsilon^(7/2) B = eval(B,epsilon^(7/2),0) -- discard epsilon^2 B = B (1 - i e epsilon / hbar phi) -- multiply by (1 - phi) B = eval(B,epsilon^(7/2),0) -- discard epsilon^2 B = B / (2 pi i hbar epsilon / m)^(3/2) -- cancel A B = B - psi() -- cancel psi B = B / epsilon -- cancel epsilon div(f) = d(f[1],x1) + d(f[2],x2) + d(f[3],x3) -- C is equation (4.18) C = hbar / i d(psi(),x) - e / c A psi() C = hbar / i div(C) - e / c dot(A,C) C = 1 / (2 m) C + e phi psi() C = -i / hbar C okay(B == C)