class CoupledFiniteDifference {

  public PDEResults1D[] solve(final CoupledPDEDataBundle pdeData1,
      final CoupledPDEDataBundle pdeData2) {
    Validate.notNull(pdeData1, "pde1 data");
    Validate.notNull(pdeData2, "pde2 data");

    final PDEGrid1D grid = pdeData1.getGrid();
    ArgumentChecker.isTrue(grid == pdeData2.getGrid(), "grids must be same object");
    final ConvectionDiffusionPDE1DCoupledCoefficients coeff1 = pdeData1.getCoefficients();
    final ConvectionDiffusionPDE1DCoupledCoefficients coeff2 = pdeData2.getCoefficients();
    final double[] initalCond1 = pdeData1.getInitialCondition();
    final double[] initalCond2 = pdeData2.getInitialCondition();
    final BoundaryCondition lower1 = pdeData1.getLowerBoundary();
    final BoundaryCondition lower2 = pdeData2.getLowerBoundary();
    final BoundaryCondition upper1 = pdeData1.getUpperBoundary();
    final BoundaryCondition upper2 = pdeData2.getUpperBoundary();
    final double lambda1 = coeff1.getLambda();
    final double lambda2 = coeff2.getLambda();

    final int tNodes = grid.getNumTimeNodes();
    final int xNodes = grid.getNumSpaceNodes();

    double[] f = new double[2 * xNodes];
    double[][] full1 = null;
    double[][] full2 = null;
    if (_showFullResults) {
      full1 = new double[tNodes][xNodes];
      full2 = new double[tNodes][xNodes];
    }
    final double[] q = new double[2 * xNodes];
    final double[][] m = new double[2 * xNodes][2 * xNodes];

    double[][] a1 = new double[2][xNodes - 2];
    final double[][] a2 = new double[2][xNodes - 2];
    double[][] b1 = new double[2][xNodes - 2];
    final double[][] b2 = new double[2][xNodes - 2];
    double[][] c1 = new double[2][xNodes - 2];
    final double[][] c2 = new double[2][xNodes - 2];

    //    final double omega = 1.5;
    //    final int oldCount = 0;
    //    final boolean omegaIncrease = false;

    double dt, t1, t2, x;
    double[] x1st, x2nd;

    System.arraycopy(initalCond1, 0, f, 0, xNodes);
    System.arraycopy(initalCond2, 0, f, xNodes, xNodes);

    if (_showFullResults) {
      if (full1 != null && full2 != null) {
        full1[0] = Arrays.copyOfRange(f, 0, xNodes);
        full2[0] = Arrays.copyOfRange(f, xNodes, 2 * xNodes);
      }
    }

    for (int i = 0; i < xNodes - 2; i++) {
      x = grid.getSpaceNode(i + 1);
      a1[0][i] = coeff1.getA(0, x);
      b1[0][i] = coeff1.getB(0, x);
      c1[0][i] = coeff1.getC(0, x);

      a1[1][i] = coeff2.getA(0, x);
      b1[1][i] = coeff2.getB(0, x);
      c1[1][i] = coeff2.getC(0, x);
    }

    final boolean first = true;
    DecompositionResult decompRes = null;

    for (int n = 1; n < tNodes; n++) {

      t1 = grid.getTimeNode(n - 1);
      t2 = grid.getTimeNode(n);
      dt = grid.getTimeStep(n - 1);

      for (int i = 1; i < xNodes - 1; i++) {

        x = grid.getSpaceNode(i);
        x1st = grid.getFirstDerivativeCoefficients(i);
        x2nd = grid.getSecondDerivativeCoefficients(i);

        q[i] = f[i];
        q[i] -= (1 - _theta) * dt * (x2nd[0] * a1[0][i - 1] + x1st[0] * b1[0][i - 1]) * f[i - 1];
        q[i] -= (1 - _theta) * dt * (x2nd[1] * a1[0][i - 1] + x1st[1] * b1[0][i - 1] + c1[0][i - 1]) * f[i];
        q[i] -= (1 - _theta) * dt * (x2nd[2] * a1[0][i - 1] + x1st[2] * b1[0][i - 1]) * f[i + 1];
        q[i] -= (1 - _theta) * dt * lambda1 * f[i + xNodes];

        q[xNodes + i] = f[xNodes + i];
        q[xNodes + i] -= (1 - _theta) * dt * (x2nd[0] * a1[1][i - 1] + x1st[0] * b1[1][i - 1]) * f[xNodes + i - 1];
        q[xNodes + i] -= (1 - _theta) * dt * (x2nd[1] * a1[1][i - 1] + x1st[1] * b1[1][i - 1] + c1[1][i - 1]) * f[xNodes + i];
        q[xNodes + i] -= (1 - _theta) * dt * (x2nd[2] * a1[1][i - 1] + x1st[2] * b1[1][i - 1]) * f[xNodes + i + 1];
        q[xNodes + i] -= (1 - _theta) * dt * lambda2 * f[i];

        a2[0][i - 1] = coeff1.getA(t2, x);
        b2[0][i - 1] = coeff1.getB(t2, x);
        c2[0][i - 1] = coeff1.getC(t2, x);

        a2[1][i - 1] = coeff2.getA(t2, x);
        b2[1][i - 1] = coeff2.getB(t2, x);
        c2[1][i - 1] = coeff2.getC(t2, x);

        m[i][i - 1] = _theta * dt * (x2nd[0] * a2[0][i - 1] + x1st[0] * b2[0][i - 1]);
        m[i][i] = 1 + _theta * dt * (x2nd[1] * a2[0][i - 1] + x1st[1] * b2[0][i - 1] + c2[0][i - 1]);
        m[i][i + 1] = _theta * dt * (x2nd[2] * a2[0][i - 1] + x1st[2] * b2[0][i - 1]);
        m[i][i + xNodes] = dt * _theta * lambda1;

        m[xNodes + i][xNodes + i - 1] = _theta * dt * (x2nd[0] * a2[1][i - 1] + x1st[0] * b2[1][i - 1]);
        m[xNodes + i][xNodes + i] = 1 + _theta * dt * (x2nd[1] * a2[1][i - 1] + x1st[1] * b2[1][i - 1] + c2[1][i - 1]);
        m[xNodes + i][xNodes + i + 1] = _theta * dt * (x2nd[2] * a2[1][i - 1] + x1st[2] * b2[1][i - 1]);
        m[xNodes + i][i] = dt * _theta * lambda2;
      }

      double[] temp = lower1.getLeftMatrixCondition(pdeData1.getCoefficients(), grid, t2);
      for (int k = 0; k < temp.length; k++) {
        m[0][k] = temp[k];
      }

      temp = upper1.getLeftMatrixCondition(pdeData1.getCoefficients(), grid, t2);
      for (int k = 0; k < temp.length; k++) {
        m[xNodes - 1][xNodes - temp.length + k] = temp[k];
      }

      temp = lower2.getLeftMatrixCondition(pdeData2.getCoefficients(), grid, t2);
      for (int k = 0; k < temp.length; k++) {
        m[xNodes][xNodes + k] = temp[k];
      }

      temp = upper2.getLeftMatrixCondition(pdeData2.getCoefficients(), grid, t2);
      for (int k = 0; k < temp.length; k++) {
        m[2 * xNodes - 1][2 * xNodes - temp.length + k] = temp[k];
      }

      temp = lower1.getRightMatrixCondition(pdeData1.getCoefficients(), grid, t1);
      double sum = 0;
      for (int k = 0; k < temp.length; k++) {
        sum += temp[k] * f[k];
      }
      q[0] = sum + lower1.getConstant(pdeData1.getCoefficients(), t2);

      temp = upper1.getRightMatrixCondition(pdeData1.getCoefficients(), grid, t1);
      sum = 0;
      for (int k = 0; k < temp.length; k++) {
        sum += temp[k] * f[xNodes - 1 - k];
      }

      q[xNodes - 1] = sum + upper1.getConstant(pdeData1.getCoefficients(), t2);

      temp = lower2.getRightMatrixCondition(pdeData2.getCoefficients(), grid, t1);
      sum = 0;
      for (int k = 0; k < temp.length; k++) {
        sum += temp[k] * f[k];
      }
      q[xNodes] = sum + lower2.getConstant(pdeData2.getCoefficients(), t2);

      temp = upper2.getRightMatrixCondition(pdeData2.getCoefficients(), grid, t1);
      sum = 0;
      for (int k = 0; k < temp.length; k++) {
        sum += temp[k] * f[xNodes - 1 - k];
      }

      q[2 * xNodes - 1] = sum + upper2.getConstant(pdeData2.getCoefficients(), t2);

      //TODO work out why SOR does not converge here
      //      final DoubleMatrix2D mM = new DoubleMatrix2D(m);
      //      final DecompositionResult res = DCOMP.evaluate(mM);
      //      f = res.solve(q);

      //      // SOR
      //
      //      int count = sor(omega, grid, freeBoundary, xNodes, f, q, m, t2);
      //      if (oldCount > 0) {
      //        if ((omegaIncrease && count > oldCount) || (!omegaIncrease && count < oldCount)) {
      //          omega = Math.max(1.0, omega * 0.9);
      //          omegaIncrease = false;
      //        } else {
      //          omega = Math.min(1.99, 1.1 * omega);
      //          omegaIncrease = true;
      //        }
      //      }
      //      oldCount = count;

      if (first) {
        final DoubleMatrix2D mM = new DoubleMatrix2D(m);
        decompRes = DCOMP.evaluate(mM);

        // first = false;
      }

      f = decompRes.solve(q);

      a1 = a2;
      b1 = b2;
      c1 = c2;

      if (_showFullResults) {
        if (full1 != null && full2 != null) {
          full1[n] = Arrays.copyOfRange(f, 0, xNodes);
          full2[n] = Arrays.copyOfRange(f, xNodes, 2 * xNodes);
        }
      }

    }
    final PDEResults1D[] res = new PDEResults1D[2];

    if (_showFullResults) {
      res[0] = new PDEFullResults1D(grid, full1);
      res[1] = new PDEFullResults1D(grid, full2);
    } else {
      final double[] res1 = Arrays.copyOfRange(f, 0, xNodes);
      final double[] res2 = Arrays.copyOfRange(f, xNodes, 2 * xNodes);
      res[0] = new PDETerminalResults1D(grid, res1);
      res[1] = new PDETerminalResults1D(grid, res2);
    }

    return res;

  }

  public  CoupledFiniteDifference();
  public  CoupledFiniteDifference(final double theta, final boolean showFullResults);

}