  Several novel imaging and non-destructive testing technologies are based on
reconstructing the spatially dependent coefficient in an elliptic partial
differential equation from measurements of its solution(s). In practical
applications, the unknown coefficient is often assumed to be piecewise constant
on a given pixel partition (corresponding to the desired resolution), and only
finitely many measurements can be made. This leads to the problem of inverting
a finite-dimensional non-linear forward operator $\mathcal F:\ \mathcal
D(\mathcal F)\subseteq \mathbb R^n\to \mathbb R^m$, where evaluating $\mathcal
F$ requires one or several PDE solutions.
  Numerical inversion methods require the implementation of this forward
operator and its Jacobian. We show how to efficiently implement both using a
standard FEM package and prove convergence of the FEM approximations against
their true-solution counterparts. We present simple example codes for Comsol
with the Matlab Livelink package, and numerically demonstrate the challenges
that arise from non-uniqueness, non-linearity and instability issues. We also
discuss monotonicity and convexity properties of the forward operator that
arise for symmetric measurement settings.
