Finite difference method: The correct formula Suppose I have a uniform 1D grid with spacing $\Delta x$ and want to solve for example the Schrödinger equation on this grid. What is the correct approximation for the second order derivative? Is it
$$
\frac{d^2\psi}{dx^2} = \frac{1}{2\Delta x^2}(\psi_{i+1} + \psi_{i-1} - 2\psi_{i})
$$ 
or is it
$$
\frac{d^2\psi}{dx^2} = \frac{1}{\Delta x^2}(\psi_{i+1} + \psi_{i-1} - 2\psi_{i})
$$ 
It would seem that I find both in literature. I personally think that the first is the correct one since it agrees with a second order Taylor expansion. On the other hand if I try to insert a plane wave of the form $\psi_{k}(x_{i}) = \text{exp}(ikx_{i})$ I obtain only the correct dispersion relation (by Taylor expansion of the cosine) for the second formula. If the last confuses you, see https://wiki.physics.udel.edu/phys824/Discretization_of_1D_Hamiltonian. 
 A: The second equation is correct. As you suggested, with
$$
\psi_{i+1} = \psi_i + \psi'_i \Delta x + \frac{1}{2} \psi''_i \Delta x^2 + \frac{1}{6} \psi'''_i \Delta x^3 + O\left(\Delta x^4\right)
$$
and
$$
\psi_{i-1} = \psi_i - \psi'_i \Delta x + \frac{1}{2} \psi''_i \Delta x^2 - \frac{1}{6} \psi'''_i \Delta x^3 + O\left(\Delta x^4\right),
$$
we have
$$
\frac{1}{\Delta x^2}\left(\psi_{i+1} + \psi_{i-1} - 2\psi_{i}\right) = \frac{\psi''_i \Delta x^2 + O\left(\Delta x^4\right)}{\Delta x^2} = \psi''_i + O\left(\Delta x^2\right)
$$
To leading order, the Hamiltonian acting on $\psi$ is then
$$
\begin{eqnarray}
\left(\hat{H} \psi\right)_i &=& -\frac{\hbar^2}{2m} \psi''_i + U_i \psi_i \\
&=& -\frac{\hbar^2}{2m \Delta x^2} \left(\psi_{i+1} + \psi_{i-1} - 2\psi_{i}\right) + U_i \psi_i \\
\end{eqnarray}
$$
As @KaneBilliot said, maybe some references have included that 2 in the denominator with the expression for $\psi''_i$?
A: Your question is essentially asking the finite different coefficients for a particular derivative order and choice of sample points. 
With this tool you can see that your second choice is correct. In general these can be found by applying a Taylor series to each term and working out coefficients that fit. 
However, in practice this is simply an algorithmic and tedious exercise and would recommend using a table or the tool I linked. 
A: You can easily derive the correct expression as follows. We can formally write down the Taylor expansion of a function as:
$$\exp(h D)f(x) = f(x+h)$$
where $D$ is the differential operator. Then we can find many possible ways to express the derivative in terms of finite differences by writing the r.h.s. in terms of finite difference operators. E.g. we can write:
$$f(x+h) = (1+\Delta)f(x)$$
where $\Delta$ acts on $f$ as:
$$\Delta f(x) = f(x+h) - f(x)$$
We can thus write:
$$\exp(h D)f(x) = (1+\Delta)f(x)$$
This allows you to formally express the differential operator in terms of the finite difference operator. 
While this will yield a formally correct expression, it will be in terms of forward differences and will thus yield an asymmetric expression. What we want here is a symmetric expression, but this obtained in just the same say, you just consider the symmetric finite difference expression:
$$\left[\exp(h D)-\exp(-h D)\right]f(x) = f(x+h) - f(x-h)$$
We can write this as:
$$\sinh(h D)f(x) = \Delta_s f(x)$$
where $\Delta_s$ is the average between the forward and backward finite difference. We can thus formally write:
$$D = \frac{1}{h}\operatorname{arcsinh}{\Delta_s} = \frac{1}{h}\left[\Delta_s - \frac{\Delta_s^3}{6 } + \frac{3 \Delta_s^5}{40 }-\frac{5 \Delta_s^7}{112 }+\cdots\right] $$
The second derivative operator is easily expressed in terms of finite differences by squaring the series, we have:
$$D^2 = \frac{1}{h^2}\left[\Delta _s^2 -\frac{\Delta _s^4}{3 }+\frac{8 \Delta _s^6}{45 } -\frac{4 \Delta _s^8}{35 }+\cdots\right]$$
Using only the first term of this series yields the desired expression:
$$\Delta_s^2 f(x) = \frac{1}{2}\Delta_s \left[f(x+h)-f(x-h)\right] = \frac{1}{4} \left[f(x+2h)-2f(x)+ f(x-2h)\right]$$
Then the smallest stepsize $h$ you can take on your grid is $h = \frac{\Delta x}{2}$, because the smallest step that appears in this formula is $2 h$. This then yields the approximation:
$$D^2 f(x)\approx \frac{f(x+\Delta x)-2f(x)+ f(x-\Delta x)}{\Delta x^2}$$
To easily find expressions for higher order terms, it's convenient to introduce the shift operator $E$ that acts like:
$$E f(x) = f(x+h) $$
Then we have:
$$\Delta_s = \frac{E - E^{-1}}{2}$$
So, the calculation of $\Delta_s^2$ given above is nothing more than the binomial expansion of the square of  the r.h.s.
These methods involving formal manipulations of differential operators and finite difference operators are also useful for purely theoretical computations, albeit they'll then yield formal expressions that will then lack a rigorous mathematical derivation. E.g. suppose that you want to compute an integral from zero to infinity of some function:
$$\int_0^{\infty} x^{s-1}f(x) dx$$
and the series expansion coefficients of the integrand are known, we have:
$$f(x) = \sum_{n=0}^{\infty}(-1)^n\frac{c_n}{n!}x^n$$
We can then write $c_n = E^{n}c_0$, where $E$ is again the shift operator. We thus have:
$$f(x) = \sum_{n=0}^{\infty}(-1)^n\frac{E^n}{n!}x^n c_0 = \exp(-E x) c_0$$
The integral can then be formally computed as:
$$\int_0^{\infty} x^{s-1}f(x) dx = \Gamma(s)E^{-s} c_0 = \Gamma(s) c_{-s}$$
Here one assumes one is allowed to analytically continue the series expansion coefficients in some way, this can be made more rigorous, see here for details. But what should be clear is that using suitably defined operators you can get to results must faster, basically by cutting through the mathematical red tape.
