I see there are some differences in the solution as I increase the resolution of my grid.
I'm using Operator Splitting to solve Diffusion Reaction equation
\begin{equation} \frac{\partial u}{\partial t} - \frac{1}{2} \nabla ^2 u = (1-u)u \end{equation}
Discretised in 1D as:
\begin{equation} (1-\frac{k}{2}A - \frac{\Delta t}{2}(1-u_E))u(t+\Delta t) = u(t) +\frac{k}{2}Au(t) + \frac{\Delta t}{2}(1-u(t))u(t) \end{equation}
where
\begin{equation} k = \frac{\Delta t}{2(\Delta x )^2} \end{equation}
\begin{equation} A= \begin{bmatrix} -2 & 1 & 0 & 0 & 0 & \ldots & 0 \\ 1 & -2 & 1 & 0 & 0 &\ldots & 0 \\ 0 & 1 & -2& 1 & 0 & \ldots & 0 \\ 0 & 0 & 1& -2 & 1 & \ldots & 0 \\ 0 & 0 & 0& 1 & -2 & \ldots & 0 \\ \vdots &\vdots &\vdots &\vdots & \ddots &\ddots &\vdots \\ 0 & 0 & 0 & \ldots &0 &1 &-2 \end{bmatrix} \end{equation}
\begin{equation} u_E= u(t) + k Au(t) + \Delta t (1-u(t))u(t) \\ \end{equation}
Is it possible to apply RK4 on an implicit scheme and improve the accuracy of the solution?
It looks to me you already invoked Forward or Backward Euler (cannot really tell, looks not correct to me) to obtain the discretized equation (there is an $\Delta t$ showing up).
What you want is actually something semi-discretized, i.e., you want to discretize the spatial derivatives only.
Here, this would imply converting the PDE
$$\partial_t u(t, x) = \frac12 \nabla^2 u(t,x) + \big(1 - u(t, x) \big) u(t,x)$$
via the second order accurate finite differences you propose to
$$U(t)' = \frac12 k A U(t) + \big(1 - U(t) \big) U(t) =: F\big(U(t)\big)$$
Now, you have just an ODE (which lacks an initial condition) for which you can use any ODE solver you like.
In particular, you could use RK4, something implicit, ...