class LocalVolPDEPricer {

  public double price(final double s0, final double k, final double r, final double b, final double t, final LocalVolatilitySurfaceStrike locVol, final boolean isCall, final boolean isAmerican,
      final int spaceNodes, final int timeNodes) {

    final double q = r - b;
    // final double s0 = fwdCurve.getSpot();
    final double sigma0 = locVol.getVolatility(t, k);
    final double mult = Math.exp(6.0 * sigma0 * Math.sqrt(t));
    final double sMin = Math.min(0.8 * k, s0 / mult);
    final double sMax = Math.max(1.25 * k, s0 * mult);

    final int tBurnNodes = (int) (USE_BURNIN ? Math.max(2, timeNodes * BURNIN_FRACTION) : 0);
    final double tBurn = USE_BURNIN ? BURNIN_FRACTION * t * t / timeNodes : 0.0;

    // set up a near-uniform mesh that includes spot and strike
    MeshingFunction xMesh = new ExponentialMeshing(sMin, sMax, spaceNodes, 0.0, new double[] {s0, k});
    MeshingFunction tMeshBurn = USE_BURNIN ? new ExponentialMeshing(0.0, tBurn, tBurnNodes, 0.0) : null;
    MeshingFunction tMesh = new ExponentialMeshing(tBurn, t, timeNodes - tBurnNodes, 0.0);
    PDEGrid1D gridBurn = USE_BURNIN ? new PDEGrid1D(tMeshBurn, xMesh) : null;
    PDEGrid1D grid = new PDEGrid1D(tMesh, xMesh);
    final int index = Arrays.binarySearch(grid.getSpaceNodes(), s0);
    ArgumentChecker.isTrue(index >= 0, "cannot find spot on grid");

    ConvectionDiffusionPDE1DStandardCoefficients coef = PDE.getBackwardsLocalVol(r, q, t, locVol);
    Function1D<Double, Double> payoff = ICP.getEuropeanPayoff(k, isCall);

    BoundaryCondition lower;
    BoundaryCondition upper;

    PDEResults1D res;

    if (isAmerican) {
      if (isCall) {
        lower = new NeumannBoundaryCondition(0.0, sMin, true);
        upper = new NeumannBoundaryCondition(1.0, sMax, false);
      } else {
        lower = new NeumannBoundaryCondition(-1.0, sMin, true);
        upper = new NeumannBoundaryCondition(0.0, sMax, false);
      }

      final Function<Double, Double> func = new Function<Double, Double>() {
        @Override
        public Double evaluate(Double... tx) {
          final double x = tx[1];
          return isCall ? Math.max(x - k, 0.0) : Math.max(k - x, 0);
        }
      };

      final FunctionalDoublesSurface free = new FunctionalDoublesSurface(func);
      if (USE_BURNIN) {
        PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> dataBurn = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(coef, payoff, lower, upper, free, gridBurn);
        PDEResults1D resBurn = INITIAL_SOLVER.solve(dataBurn);

        PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> data = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(coef, resBurn.getTerminalResults(), lower, upper, free, grid);
        res = SOLVER.solve(data);
      } else {
        PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> data = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(coef, payoff, lower, upper, free, grid);
        res = SOLVER.solve(data);
      }

    } else {
      if (isCall) {
        lower = new NeumannBoundaryCondition(0.0, sMin, true);
        Function1D<Double, Double> upFunc = new Function1D<Double, Double>() {
          @Override
          public Double evaluate(Double time) {
            return Math.exp(-q * time);
          }
        };
        upper = new NeumannBoundaryCondition(upFunc, sMax, false);
      } else {
        Function1D<Double, Double> downFunc = new Function1D<Double, Double>() {
          @Override
          public Double evaluate(Double time) {
            return -Math.exp(-q * time);
          }
        };
        lower = new NeumannBoundaryCondition(downFunc, sMin, true);
        upper = new NeumannBoundaryCondition(0.0, sMax, false);
      }
      if (USE_BURNIN) {
        PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> dataBurn = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(coef, payoff, lower, upper, gridBurn);
        PDEResults1D resBurn = INITIAL_SOLVER.solve(dataBurn);

        PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> data = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(coef, resBurn.getTerminalResults(), lower, upper, grid);
        res = SOLVER.solve(data);
      } else {
        PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients> data = new PDE1DDataBundle<ConvectionDiffusionPDE1DCoefficients>(coef, payoff, lower, upper, grid);
        res = SOLVER.solve(data);
      }
    }

    return res.getFunctionValue(index);
  }

}