I am implementing a finite difference scheme for the following PDE: $$ u_t=u_{xx}+f(u) $$ with the nonflux boundary condition and a given initial condition. The scheme is backward in time and central in space(which should be unconditionally stable) as $$ \frac{u^{k+1}_i-u^k_i}{\Delta t}=\frac{u^{k+1}_{i-1}-2u^{k+1}_{i}+u^{k+1}_{i+1}}{\Delta x^2}+f(u_i^k) $$ however, I am not getting the correct order of convergence, i.e. first order in time and second order in space. And the scheme is not convergence(actually blows up...). I am just wondering is there any special condition I need to be careful on reaction-diffusion equations? because I downgrade the case into heat equation $$ u_t=u_{xx} $$ with the same boundary condition and then I am getting the expected order of convergence.
The update scheme is $$ A\mathbf{u}^{k+1} = u^{k+1}_i - \alpha (u^{k+1}_{i-1}-2u^{k+1}_{i}+u^{k+1}_{i+1}) = u^k_i + \Delta t f(u_i^k) $$ where $A$ has diagonal = $1+2\alpha$ with off diagonals = $-alpha$ and the $A(1,2)=A(N,N-A)=-2\alpha$. For the heat equation, there's no $\Delta t f(u_i^k)$ on the right-hand side.
To take care of the nonflux boundary condition, I also used second-order central in space as $$ \frac{u_{-1}-u_1}{2\Delta x}=0 $$
I do suspect the boundary is somehow wrong because I plotted each time stamp to see how spatial solution evolves and I do observe errors getting larger at the right end. But I can't think of other ways to take care of the boundaries... any suggestion is appreciated!
update: please please update MATLAB version. I used the same scheme but with the latest MATLAB version and got the right order of convergence...