I modeled the heat equation,$$u_t = au_{xx}$$using the common 2nd-order Crank-Nicolson scheme,$$\frac{u^{n+1}_i-u^{n}_i}{dt} = \frac{a}{2\,dx}\left(u_{i-1}^{n+1}+u_{i+1}^{n+1}-2u_i^{n+1} + u_{i-1}^n+u_{i+1}^n-2u_i^n\right)$$and the standard 4th-order Crank-Nicolson scheme,$$\frac{u^{n+1}_i-u^{n}_i}{dt} = \frac{a}{24\,dx}\left(-u_{i-2}^{n+1}+16u_{i-1}^{n+1}-30u_{i}^{n+1}+16u_{i+1}^{n+1}-u_{i+2}^{n+1}\right) \\   \quad+\frac{a}{24\,dx}\left(-u_{i-2}^{n}+16u_{i-1}^{n}-30u_{i}^{n}+16u_{i+1}^{n}-u_{i+2}^{n}\right)$$Outside of figuring out how to solve the pentadiagonal matrix (without resorting to libraries like LAPACK), this was fairly straight-forward.But my results were strange: the 2nd-order scheme was more accurate to the exact solution (in 3 different test cases!) than the 4th-order scheme. Is there any reason we should expect a 2nd-order accurate scheme to perform better than the 4th-order scheme?If it matters, $a=1$, $x\in(0,1)$ and $u(0,t)=u(1,t)=0$ for all examples:$u(x,0)=6\sin(\pi x)$ leading to the analytic solution $u(x,t)=6\sin(\pi x)\exp(-\pi^2t)$$u(x,0)=12\sin(9\pi x) - 7\sin(4\pi x)$ leading to $u(x,t)=12\sin(9\pi x)\exp(-9\pi^2t) - 7\sin(4\pi x)\exp(-9\pi^2t)$$u(x,0)=10$ leading to $u(x,t)=\sum_{k=1}^{100}\frac{40}{(2k-1)\pi}\exp\left(-(2k-1)^2\pi^2t\right)\sin\left((2k-1)\pi x\right)$With this last one, though, both did pretty poorly at the boundaries, but the 2nd-order method was closer to the exact in the range $x\in(0.1,0.9)$.EDITYou can see that it's an order of magnitude difference between the two schemes.two sine termsingle sine termconstant 10EDIT 2I ran this with different $N$'s and still found the 2nd order version to be more accurate than the 4th order. 
Crank-Nicolson for 2nd- and 4th-order finite differences
finite difference;crank nicolson