final PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> db2 = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(pde2, initalCondition, lower, upper, grid);

    final ThetaMethodFiniteDifference thetaMethod = new ThetaMethodFiniteDifference(0.5, true);

    final PDEFullResults1D res2 = (PDEFullResults1D) thetaMethod.solve(db2);

    PDEUtilityTools.printSurface("State 1 density", res1);

    PDEUtilityTools.printSurface("State 2 density", res2);

  }

    final ThetaMethodFiniteDifference thetaMethod = new ThetaMethodFiniteDifference(0.5, true);

    final PDEFullResults1D res1 = (PDEFullResults1D) thetaMethod.solve(db1);

    PDEUtilityTools.printSurface("State 1 density", res1);

    PDEUtilityTools.printSurface("State 2 density", res2);

  }

    final MeshingFunction spaceMesh = new HyperbolicMeshing(minMoneyness, maxMoneyness, centreMoneyness, nStrikeSteps, strikeMeshBunching);

    final PDEGrid1D grid = new PDEGrid1D(timeMesh, spaceMesh);

    final Function1D<Double, Double> intCond = (new InitialConditionsProvider()).getForwardCallPut(isCall);

    return res;

  }

  private PDEFullResults1D runForwardPDESolver(final ForwardCurve forwardCurve, final LocalVolatilitySurfaceStrike localVolatility,

      final boolean isCall, final double theta, final double maxT, final double maxAbsProxyDelta, final int nTimeSteps, final int nStrikeSteps,

    //keep the grid the same regardless of spot (useful for finite-difference)

    final MeshingFunction spaceMesh = new HyperbolicMeshing(0.0, maxFwd, fwdNodeCentre, nFwdNodes, spotMeshBunching);

    final PDEGrid1D grid = new PDEGrid1D(timeMesh, spaceMesh);

    final PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> db = new PDE1DDataBundle<>(pde, payoff, lower, upper, grid);

    return res;

  }

  /**

   * Convert the results of running the forward PDE, which are forward option prices divided by the relevant forward, to an implied volatility

    final BoundaryCondition upper = new NeumannBoundaryCondition(1.0, upperLevel, false);

    final MeshingFunction timeMesh = new ExponentialMeshing(0.0, EXPIRY, nTimeNodes, 6.0);

    final MeshingFunction spaceMesh = new HyperbolicMeshing(0, upperLevel, STRIKE, nSpotNodes, 0.05);

    final PDEGrid1D grid = new PDEGrid1D(timeMesh, spaceMesh);

    final int fwdIndex = grid.getLowerBoundIndexForSpace(forward);

    final double[] fwd = new double[4];

    final double[] vol = new double[4];

    for (int i = 0; i < 4; i++) {

    final MeshingFunction timeMeshF = new ExponentialMeshing(0, t, 50, 4.0);

    final MeshingFunction timeMeshB = new DoubleExponentialMeshing(0, t, t / 2, 50, 2.0, -4.0);

    final PDEGrid1D grid = new PDEGrid1D(timeMeshF, spaceMeshF);

    PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> dbF = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(pde, initialCond, lower, upper, grid);

    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);

      BoundaryCondition upperB = new NeumannBoundaryCondition(isCall ? 1 : 0, sH, false);

      Function1D<Double, Double> bkdIC = initialConProvider.getEuropeanPayoff(k, isCall);

      PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> dbB = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(pdeB, bkdIC, lowerB, upperB, gridB);

      double price1 = resB.getFunctionValue(index);

      double price2 = resB.getFunctionValue(index + 1);

      double vol1 = BlackFormulaRepository.impliedVolatility(price1, s1, k, t, isCall);

      double vol2 = BlackFormulaRepository.impliedVolatility(price2, s2, k, t, isCall);

      double volBPDE = w * vol1 + (1 - w) * vol2;

 

    //run the PDE solver backward to the dividend date

    // PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> db1 = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(pde, initialCon, lower1, upper1, grid1);

    final PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> db1 = new PDE1DDataBundle<>(pde, payoff, lower, upper, grid);

    final Interpolator1DDataBundle interpolDB = INTEPOLATOR1D.getDataBundle(sNodes, res.getTerminalResults());

    final double val = INTEPOLATOR1D.interpolate(interpolDB, lnFT);

    assertEquals(0.41491529, Math.sqrt(-2 * (val) / EXPIRY), 5e-4); //Number from backwardsPDETest

    final MeshingFunction spaceMesh = new ExponentialMeshing(yL, yH, 101, 0.0);

    final PDEGrid1D grid = new PDEGrid1D(timeMesh, spaceMesh);

    final PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> db = new PDE1DDataBundle<>(pde, PURE_LOG_PAY_OFF, lower, upper, grid);

    final int n = res.getNumberSpaceNodes();

    final double val = res.getFunctionValue(n / 2);

    return val;

    final double[] sNodes1 = grid1.getSpaceNodes();

    //run the PDE solver backward to the dividend date

    final PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> db1 = new PDE1DDataBundle<>(pde, initialCon, lower1, upper1, grid1);

    //Map the spot nodes after (in calendar time) the dividend payment to nodes before

    final int nSNodes = sNodes1.length;

    final double[] sNodes2 = new double[nSNodes];

    final double lnBeta = Math.log(1 - BETA);