Noninvasive imaging based on wave scattering remains a difficult problem in those cases where the forward map can only be adequately simulated by solving the appropriate partialdifferential equation subject to boundary conditions. We develop a method for solving these linear boundary-value problems which is efficient and exact, trading off storage requirements against computation time. The method is based on using the present solution within the Woodbury formula for updating solutions given changes in the trial image, or state. Hence the method merges well with the Metropolis-Hastings algorithm using localized updates. The scaling of the method as a function of image size and measurement set size is given. We conclude that this method is considerably more efficient than earlier algorithms that we have used to demonstrate sampling for inverse problems in this class. We give examples of sampling for imaging electrical conductivity from a simple synthetic data set. Full Bayesian inference is demonstrated with expectations calculated over the posterior for Potts-type prior distributions.