A few years ago Andreasen and Huge have introduced an efficient and arbitrage free volatility interpolation method  based on a one step finite difference implicit Euler scheme applied to a local volatility parametrization. Probably the most notable use case is the generation of a local volatility surface from a set of option quotes.
Starting point is Dupire’s forward equation for the prices of European call options at time with strike and maturity is given by
Define the normalized call price in terms of the discount factor , the forward price and the moneyness as
The Dupire forward equation for the normalized prices is then given by
Rewriting this equation in terms of amd yields
The normalized put prices are fulfilling the same equation, which can easily been shown by inserting the call-put parity into the equation above
The numerical stability of the original algorithm  can be enhanced for deep ITM options by calibrating to calls and puts instantaneously. Also the interpolation scheme has a significant impact on the stability. This topic has been discussed in . Using concentrated meshes along the current spot level for the finite difference scheme is of advantage for the stability and accuracy of the algorithm.
In order to stabilize the calculation of the local volatility function
one should evaluate the first order derivative of w.r.t. time using the fact that the derivative of the inverse of the matrix is given by
As an example, the diagram below shows different calibrations of the Andreasen-Huge volatility interpolation to a SABR volatility skew at discrete strike sets
for the SABR parameter
The QuantLib implementation is part of the pull request #380.
 J. Andreasen, B. Huge, Volatility Interpolation
 F. Le Floc’h, Andreasen-Huge interpolation – Don’t stay flat
 J. Healy, A spline to fill the gaps with Andreasen-Huge one-step method