function [x, y, u, f] = iteration(iterations, N)

  dphi = pi / (2 * (N - 1));
  phi = (0:N - 1) * dphi;
  N = columns(phi);
  x = sin(phi);
  x = repmat(x, N, 1);
  y = sqrt(0.5) * sin(phi)';
  y = repmat(y, 1, N);
  
  dx = dphi * cos(phi);
  dx = repmat(dx, N, 1);
  dy = dphi * sqrt(0.5) * cos(phi)';
  dy = repmat(dy, 1, N);
  
  u = y;
    
  for i = 1:iterations    
    jacobi = jacobi_function(u, dx, dy);
    f = f_function(u, dx, dy);
    
    du = -f ./ jacobi;
    
    du(logical(flipud(eye(N)))) = 0;
    
    u = u + du;
  endfor
  
  u(logical(flipud(triu(ones(N, N)) - eye(N)))) = 0;
  f = f_function(u, dx, dy);
  f(logical(flipud(triu(ones(N, N)) + diag(ones(1, N - 1), -1)))) = 0;
  
endfunction
