I was trying to solve numerically the classic Black-Scholes (B-S) problem $$\frac{\partial{u}}{\partial{\tau}}=\frac{1}{2}\sigma^2x^2\frac{\partial^2{u}}{\partial{x}^2}+rx \frac{\partial{u}}{\partial{x}}-ru,~x\in (0,\infty),~\tau\in(0,T),$$ with $u(x,0)=f(x)$ and $u(0,\tau)=g(\tau),~~u(x,\tau)\sim h(\tau)$ as $x$ tends to $\infty$. The Crank-Nicoson scheme was used for time integration of system of ODE $$\begin{array}{llll} \displaystyle{\frac{\partial{u_i}(\tau)}{\partial{\tau}}}&=&\displaystyle{\frac{1}{2h^2}\sigma^2x_i^2(u_{i+1}(\tau)-2u_i(\tau)+u_{i-1}(\tau))}& \\ & & \displaystyle{+\frac{1}{2h}rx_i(u_{i+1}(\tau)-u_{i-1}(\tau))-ru_i(\tau)},& \\ & = & \displaystyle{(\frac{1}{2h^2}\sigma^2x_i^2+\frac{1}{2h}rx_i)u_{i+1}(\tau)+(-\frac{1}{h^2}\sigma^2x_i^2-r)u_i(\tau)}& \\ & &\displaystyle{+(\frac{1}{2h^2}\sigma^2x_i^2-\frac{1}{2h}rx_i)u_{i-1}(\tau)},~~i=2,3,...,N-1,& \end{array} $$
formed by the space discretization using classical finite differences. Interestingly the error data found as

where the symbols are described by: Let $\mathcal U_{i,j}^{M,N}$ be the numerical approximation with $M$ and $N$ points in space and time directions respectively, The closed form of exact solution for the above B-S problem is given by \begin{equation}\label{eq23} u(x, t)= \exp (x) \mathcal N\left(\hat{d}_{1}\right)-E \exp (-\hat{r} t) \mathcal N\left(\hat{d}_{2}\right), \end{equation} where $$ \hat{d}_{1}=\frac{x-\ln E+\left(\hat{r}+\frac{1}{2} \hat{\sigma}^{2}\right) t}{\hat{\sigma} \sqrt{t}} ~~\text{and}~~\hat{d}_{2}=\hat{d}_{1}-\hat{\sigma} \sqrt{t}.$$ Compute the true $L_{\infty}$ norm error (maximum absolute error, $e_{\max}^{M,N}$), $L_{2}$ norm error (root mean square error, $e_{rms }^{M,N}$) and corresponding orders of convergence $p_{max }^{M,N}$ and $p_{rms }^{M,N}$ as follows: $$e_{\max }^{M, N}=\max _{0 \leq m \leq M}\left|u\left(x_{m}, t_{N}\right)-\mathcal U_{m,N}^{M,N}\right|,$$ $$e_{rms}^{{M,N}}=\sqrt{\frac{\sum_{m=0}^{M}\left(u\left(x_{m}, t_{N}\right)-\mathcal U_{m,N}^{M,N}\right)^{2}}{M+1}},$$ and $$p_{\max}^{M,N}=\log_2\left(\frac{e_{\max}^{M,N}}{e_{\max}^{2M,2N}}\right),$$ $$p_{rms}^{M,N}=\log_2\left(\frac{e_{rms}^{M,N}}{e_{rms}^{2M,2N}}\right).$$ It is observed from the data that the convergence order of Crank-Nicolson method is only 1.0. How can we explain the order comes like apart from 2?