Метод релаксации неправильно работает

void Relax(int N, double Omega, double eps, double phi[LEN][LEN], PF f, double U[LEN][LEN])
{
 double h = 3/(double)N;
 double ro[LEN][LEN];
 double c[LEN][LEN];
 for (int i = 0; i < N; i++)
  {
   double x = i * h;
   for (int j = 0; j < N; j++)
    {
      double y = j * h;
      if (x >= 1 && x <= 2)
       {
        if (y >= 1 && y <= 2)
         ro[i][j] = f(x,y);
        else
         ro[i][j] = 0;
       }
      else ro[i][j] = f(x,y);
    }
  }
 
  for (int i = 1; i < N-1; i++)
    for (int j = 1; j < N-1; j++)
     {
      if (i * h >= 1 && i * h <= 2)
       if (j*h >= 1 && j*h <= 2)
        U[i][j] = (1-Omega) * phi[i][j] + (Omega * (phi[i+1][j] + phi[i-1][j] + phi[i][j+1] + phi[i][j-1] + h*h*ro[i][j]))/4;
       else U[i][j] = 0;
      else  U[i][j] = (1-Omega) * phi[i][j] + (Omega/4) * (phi[i+1][j] + phi[i-1][j] + phi[i][j+1] + phi[i][j-1] + h*h*ro[i][j]);
     }
  Diff(U,phi,c,N);
 
  int k = 0;
 
  while (Norm1(c,N) >= eps)
  {
   for (int i = 0; i < N; i++)
    for (int j = 0; j < N; j++)
     phi[i][j] = U[i][j];
 
  for (int i = 1; i < N-1; i++)
    for (int j = 1; j < N-1; j++)
     {
      if (i * h >= 1 && i * h <= 2)
       if (j*h >= 1 && j*h <= 2)
        U[i][j] = (1-Omega) * phi[i][j] + (Omega/4) * (phi[i+1][j] + phi[i-1][j] + phi[i][j+1] + phi[i][j-1] + h*h*ro[i][j]);
       else U[i][j] = 0;
      else  U[i][j] = (1-Omega) * phi[i][j] + (Omega/4) * (phi[i+1][j] + phi[i-1][j] + phi[i][j+1] + phi[i][j-1] + h*h*ro[i][j]);
     }
 
   Diff(U,phi,c,N);
   Form1->ProgressBar1->Position = ++k;
  }
}