final BoundaryCondition lower = new DirichletBoundaryCondition(1.0, xL);

    final BoundaryCondition upper = new NeumannBoundaryCondition(0.0, xH, false);

    final MeshingFunction spaceMesh = new HyperbolicMeshing(xL, xH, 1.0, 40, 0.05);

    final MeshingFunction timeMesh = new ExponentialMeshing(0, 2.0, 30, 0.2);

    final Function1D<Double, Double> initialCond = INITIAL_COND_PROVIDER.getForwardCallPut(true);

    final Function1D<DoubleMatrix1D, DoubleMatrix1D> volFunc = new Function1D<DoubleMatrix1D, DoubleMatrix1D>() {

      @Override

    final BoundaryCondition lower = new DirichletBoundaryCondition(0, 0);

    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 PDEResults1D res = solver.solve(new PDE1DDataBundle<>(pde, payoff, lower, upper, grid));

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

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

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

      final double price = res.getFunctionValue(i + fwdIndex - 1);

      vol[i] = BlackFormulaRepository.impliedVolatility(price, fwd[i], STRIKE, EXPIRY, true);

    }

    final Interpolator1DDoubleQuadraticDataBundle idb = INTERPOLATOR_1D.getDataBundle(fwd, vol);

    final int tNodes = 51;

    final int xNodes = 151;

    final MeshingFunction timeMesh = new ExponentialMeshing(0, T, tNodes, 7.5);

    final MeshingFunction spaceMesh = new HyperbolicMeshing(0, 10 * SPOT, SPOT, xNodes, 0.01);

    final PDEFullResults1D res = PRICER.solve(grid, theta);

    final double[] expiries = timeMesh.getPoints();

    final double[] strikes = spaceMesh.getPoints();

    final double[] forwards = new double[tNodes];

    BoundaryCondition upper = new NeumannBoundaryCondition(0.0, xH, false);

    final MeshingFunction spaceMeshF = new HyperbolicMeshing(xL, xH, 1.0, 200, 0.001);

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

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

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

    final MeshingFunction timeMesh = new ExponentialMeshing(0, EXPIRY, 100, 0.0);

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

    //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 PDETerminalResults1D res = (PDETerminalResults1D) solver.solve(db1);

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

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

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

    final PDEResults1D res = solver.solve(db);

    final int n = res.getNumberSpaceNodes();

    final MeshingFunction timeMesh1 = new ExponentialMeshing(0, EXPIRY - DIVIDEND_DATE - 1e-6, 50, 0.0);

    final MeshingFunction timeMesh2 = new ExponentialMeshing(EXPIRY - DIVIDEND_DATE + 1e-6, EXPIRY, 50, 0.0);

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

    //run the PDE solver backward to the dividend date

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

    final PDETerminalResults1D res1 = (PDETerminalResults1D) solver.solve(db1);

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

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

      final double temp = sNodes1[i];

      if (temp < 0) {

        sNodes2[i] = Math.log(Math.exp(temp) + ALPHA) - lnBeta;

      }

      else {

        sNodes2[i] = temp + Math.log(1 + ALPHA * Math.exp(-temp)) - lnBeta;

      }

    }

    final BoundaryCondition lower2 = new NeumannBoundaryCondition(1.0, sNodes2[0], true);

    final BoundaryCondition upper2 = new NeumannBoundaryCondition(1.0, sNodes2[nSNodes - 1], false);

    //run the PDE solver backward from the dividend date to zero

    final PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> db2 = new PDE1DDataBundle<>(pde, res1.getTerminalResults(), lower2, upper2, grid2);

    final int tNodes = 20;

    final int xNodes = 100;

    final MeshingFunction timeMesh = new ExponentialMeshing(0, 5, tNodes, 5.0);

    final MeshingFunction spaceMesh = new HyperbolicMeshing(0, 6 * SPOT, SPOT, xNodes, 0.01);

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

      mc.solve(grid, 0.5);

    }

    if (benchmarkCycles > 0) {

    final int xNodes = 100;

    final MeshingFunction timeMesh = new ExponentialMeshing(0, t, tNodes, 2.0);

    final MeshingFunction spaceMesh = new HyperbolicMeshing(0, 6.0 * FORWARD_CURVE.getForward(t), FORWARD_CURVE.getSpot(), xNodes, 0.01);

    //TwoStateMarkovChainDensity densityCal = new TwoStateMarkovChainDensity(forward, chainData);

    final TwoStateMarkovChainWithLocalVolDensity densityCal = new TwoStateMarkovChainWithLocalVolDensity(FORWARD_CURVE, DATA, new AbsoluteLocalVolatilitySurface(ConstantDoublesSurface.from(1.0)));

    final PDEFullResults1D[] denRes = densityCal.solve(grid);

    System.out.println("Densities ");

    PDEUtilityTools.printSurface("state 1 density", denRes[0]);

    MeshingFunction timeMesh = new ExponentialMeshing(0, T, tNodes, 5.0);

    //MeshingFunction spaceMesh = new ExponentialMeshing(0.0, 6.0 * FORWARD.getForward(T), xNodes, 3.0);

    MeshingFunction spaceMesh = new HyperbolicMeshing(-0.0 * FORWARD.getForward(T), 6.0 * FORWARD.getForward(T), FORWARD.getSpot(), xNodes, 0.01);

    PDEFullResults1D[] res = DENSITY_CAL.solve(grid);

    PDEUtilityTools.printSurface("state 1 density", res[0]);

    PDEUtilityTools.printSurface("state 2 density", res[1]);