numerical solution should be improved with the mesh refinement, but it didn't. How should I interpret those results?
Unfortunately, all one can say is it means something's wrong.
Reducing $\Delta x$ should decrease the error (compared to the analytic solution). This is because techniques for numerically estimating derivatives get more accurate with smaller $\Delta x$. If your code was perfectly correct, then the only error would be this numerical error (i.e. error due to the grid and timestep not being infinitesimally fine). But because your result isn't improving with finer grids in this way, it shows there's "some other problem".
You can analytically verify that the given analytical solution really is a solution to the governing differential equations (i.e. substitute the solution into the governing equation, and work out the derivatives, and then you can confirm that you end up with those given source terms). A symbolic calculator is probably needed, like
sympy (for python).
So, there's something wrong with your program - a bug. It could be anything, from the cfd code (e.g. how you discretize), to typos in your program.
Apart from standard debugging, code review etc, one way to narrow it down is to simplify the governing equation, by removing some terms. Then, substitute the same given analytic solution from before into the governing equation - which should yield different source terms. (Because the solution is the same as before, the initial and boundary conditions stay the same). Now, do as you did before, but with the new source terms: run the simulation, and compare the result with the analytic solution. Each time you refine the grid-size (i.e. reduce $\Delta x$), the error should decrease. By seeing which of the simplified governing equations behave coreectly or incorrectly, you can narrow down which terms are causing the problem, and thus which part of the code and/or program are buggy.
BTW This is the Method of Manufactured Solutions, though there's usually more care with how much the error decreases with decreasing $\Delta x$ - it should be $O( (\Delta x)^p )$, where $p$ is the order of the estimation technique (for example finite difference "forward Euler" is 1st order, $p=1$; central difference is 2nd order, $p=2$).
Because programs usually have bugs, it's unlikely that this behaviour is not caused by a bug. But, if the code and your program are 100% bug-free, then it could be a deeper issue with this scheme and/or this solver, for this particular problem. But it's almost certainly a bug. So look into that first.