There are several possible reasons for the solution to not match the analytical result that I could see just by looking at the code.
This seems to me a forward-euler (first order explicit) discretization of the ODE. This may not be the best choice and in-fact may not even be stable for larger more complicated problems. Try implementing backward euler for instance.
Also there is a single-pass where you solve for Ci and then update u. If you are solving two coupled equations in a non-monolithic fashion then you cannot do this. See for instance Jeremy Bleyer’s tutorial
This in addition to my previous comment on using an inf-sup stable mixed element, such as Taylor-Hood (\mathbb{P}_k-P_{k-1}) element should be the way to go.