Approximating Schrodinger's equation
Recall that:
|Ψ(t)>=e^(-iHT/h) |Ψ(0)>
A better approximation is to use the Cayley form.
e^(-iHt/h) ~= (1- (1/2) i t H)/(1+ (1/2) i t H)


Now we can use the finite-difference equation we used before to come up with:
This is a better equation to find the wave function at a later time, but
it is implicit rather than an explicit formula.
(Note that this method assumes the potential is constant in time).
Aside: What is going on here?
One way to look at this is in matrix form:
|Ψ(t+dt)>=U(dt)|Ψ(t)>
In this picture, Ψ is a vector and U is a unitary matrix.
(Note that this is similar to a Markov process).
In this first-order Cayley approximation, we know that the wave function
at t depends on only the adjacent values of x. That is, the matrix is in
the following "tridiagonal" form:

Note that this is just a system of coupled equations.
Using either the matrix above or the equation developed earlier, we can
turn this into a system of equations. Since the matrix is tridaigonal,
each equation will have three terms (Ψ at x-1,x,x+1). If we assume
that &Psi at the endpoints (-L, L) is constrained, we can iteratively
solve the equation. After lots of mathematical chug and plug we can
get:
Ψ(L,t+dt) = -f[L,t] / e[L]
Ψ(x,t+dt) = (Ψ(x+dx,t+dt)-f[x,t]) / e[x-dx]
(e and f are terms that are functions of the potential and the Ψ at
the previous time step)
That is, we get a recurrence so that, if we iterate from x=L down to
x=-L, we can get all values of the wave function.
Results!