from what I can tell, you want to monitor the value of g as a function of mesh resolution, but it cannot be plotted because I think it is not indexed correctly. can you show what your [non-]convergence looks like? It's possible since you appear to be using first order derivatives that you just aren't getting high enough resolution to get into the convergence regime.
how are you assessing convergence? are you trying to fit a power law in mesh size?
also your code is inefficient - pay attention to teh warnings about preallocation and sparse indexing.
