PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> dbF = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(pde, initialCond, lower, upper, grid);
PDETerminalResults1D res = (PDETerminalResults1D) solver.solve(dbF);
final double minK = Math.exp(-6 * rootT);
final double maxK = Math.exp(6 * rootT);
Map<Double, Double> vols = PDEUtilityTools.priceToImpliedVol(fwdCurve, t, res, minK, maxK, true);
DoubleQuadraticInterpolator1D interpolator = Interpolator1DFactory.DOUBLE_QUADRATIC_INSTANCE;
Interpolator1DDataBundle idb = interpolator.getDataBundle(vols);
//set up for solving backwards PDE
ConvectionDiffusionPDE1DStandardCoefficients pdeB = pdeProvider.getBackwardsLocalVol(t, lvsm);
double sL = xL * spot;
double sH = xH * spot;
final MeshingFunction spaceMeshB = new HyperbolicMeshing(sL, sH, spot, 200, 0.001);
final PDEGrid1D gridB = new PDEGrid1D(timeMeshB, spaceMeshB);
int index = SurfaceArrayUtils.getLowerBoundIndex(gridB.getSpaceNodes(), spot);
double s1 = gridB.getSpaceNode(index);
double s2 = gridB.getSpaceNode(index + 1);
final double w = (s2 - spot) / (s2 - s1);
//solve a separate backwards PDE for each strike
for (int i = 0; i < 10; i++) {
double z = -5 + i;
double k = spot * Math.exp(0.4 * rootT * z);
double x = k / fwd;
double vol = ivs.getVolatility(t, k);
double volFPDE = interpolator.interpolate(idb, x);
boolean isCall = (k >= fwd);
BoundaryCondition lowerB = new NeumannBoundaryCondition(isCall ? 0 : -1, sL, true);
BoundaryCondition upperB = new NeumannBoundaryCondition(isCall ? 1 : 0, sH, false);