The “solution” of an inverse problem may be thought of as the parameter choice that “best” fits the observations. Allowing for uncertainty in the observations and sensitivity of the data-to-parameter map suggest that it’s appropriate to instead think of the solution of the inverse problem in terms of (posterior) probability distribution for the sought parameter, such as provided by a Bayesian perspective [1]. A formal solution for the a posteriori probability distribution (PPD) for inverse problems is often readily available, but is difficult to work with practically. For unimodal posteriors, the PPD may be approximately characterized by its mode and covariance about the mode. Computing these quantities can present formidable challenges. In practice, finding the mode may require minimizing a very high dimensional function (O(10^5−10^6) parameters) with complex (e.g. discretized nonlinear) PDE constraints. For such a problem, merely storing the final covariance may require tens of terabytes. Nevertheless, it is known that the structure of inverse problems leads to a sparse update on a predictable and a priori known covariance operator. We propose a modification of a popular Newton-CG algorithm to find the mode point of the PPD. This simple adaptation takes advantage of the structure of the inverse problem and uses results from computational linear algebra to compute the covariance while finding the mode point. In particular, computing the covariance requires no additional PDE solves, and thus its computational cost is relatively “free.” We present the algorithm, prove a theorem that shows that the algorithm converges in exact arithmetic, and present several numerical examples that show its performance. Computations are performed using the open source libraries hIPPYlib [2] and FEniCS [3].