  We present a numerical scheme that can be combined with any fixed boundary
finite element based Poisson or Grad-Shafranov solver to compute the first and
second partial derivatives of the solution to these equations with the same
order of convergence as the solution itself. At the heart of our scheme is an
efficient and accurate computation of the Dirichlet to Neumann map through the
evaluation of a singular volume integral and the solution to a Fredholm
integral equation of the second kind. Our numerical method is particularly
useful for magnetic confinement fusion simulations, since it allows the
evaluation of quantities such as the magnetic field, the parallel current
density and the magnetic curvature with much higher accuracy than has been
previously feasible on the affordable coarse grids that are usually
implemented.
