I suggest you try the following:

C=A+iB, solve for Ax=b and Bx=d, where RHS y=b+id, then try to look for intersection of solution sets associated with the two above seperate real linear systems.

solve for Mx=Y, where M is 2m by n matrix, Y is 2m by 1 vector such that M =[ A , B ], where first m rows of M are the rows of A (m by n) and next m rows are the rows of B (m by n) , likewise Y =[ b , d] where the first m rows are the b vector and next m rows are the d vector. a solution of the above linear system, if exists, is a real valued solution to the original linear system Cx=y.

You can do a test. You obviously do not know the solution to the original problem, but propose a solution, substitute it into the PDE, and boundary conditions and add appropriate source terms to both the PDE and boundary conditions such that your proposed solution satisfies the new PDE and boundary conditions. After this, you can do a convergence study (well behaving sources do not effect convergence orders). I haven't read the paper, but based on the elements you use, and what function space the solution exists in, you can figure out the convergence order as you refine the mesh. I think in your case, you would find that you are not yielding the correct convergence order (i.e there is a bug in your code).

Thanks dear Ryan Vogt and Debopam Ghosh