function jacobi = jacobi_function(u, dx, dy)
  
  [u_plus_null, u_mins_null, u_null_plus, u_null_mins] = shift_values(u);
  
  # derivaltak
  u_x = (u_plus_null - u_mins_null) ./ (2 * dx);
  u_y = (u_null_plus - u_null_mins) ./ (2 * dy);
  u_x_x = (u_plus_null + u_mins_null - 2 * u) ./ (dx .^ 2);
  u_y_y = (u_null_plus + u_null_mins - 2 * u) ./ (dy .^ 2);
  
  [u_x_plus_null, u_x_mins_null, u_x_null_plus, u_x_null_mins] = shift_values(u_x);
  [u_y_plus_null, u_y_mins_null, u_y_null_plus, u_y_null_mins] = shift_values(u_y);
  
  # g
  g = (1 + u_x .^ 2 + u_y .^ 2) .^ -0.5;
  gd_plus_null = (1 + u_x_plus_null .^ 2 + u_y_plus_null .^ 2) .^ -1.5;
  gd_mins_null = (1 + u_x_mins_null .^ 2 + u_y_mins_null .^ 2) .^ -1.5;
  gd_null_plus = (1 + u_x_null_plus .^ 2 + u_y_null_plus .^ 2) .^ -1.5;
  gd_null_mins = (1 + u_x_null_mins .^ 2 + u_y_null_mins .^ 2) .^ -1.5;

  # fuggveny visszateresi ertek
  
  jacobi = 0.25 * (gd_plus_null + gd_mins_null) .* u_x .^ 2 + 0.25 * (gd_null_plus + gd_null_mins) .* u_y .^ 2 - 2 * g .* (1 ./ dx .^ 2 + 1 ./ dy .^ 2);
  
endfunction