function f = f_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_plus_null_plus_null, u_plus_null_mins_null, u_plus_null_null_plus, u_plus_null_null_mins] = shift_values(u_plus_null);
  [u_mins_null_plus_null, u_mins_null_mins_null, u_mins_null_null_plus, u_mins_null_null_mins] = shift_values(u_mins_null);
  [u_null_plus_plus_null, u_null_plus_mins_null, u_null_plus_null_plus, u_null_plus_null_mins] = shift_values(u_null_plus);
  [u_null_mins_plus_null, u_null_mins_mins_null, u_null_mins_null_plus, u_null_mins_null_mins] = shift_values(u_null_mins);
  
  u_x_plus_null = (u_plus_null_plus_null - u_mins_null_plus_null) ./ (2 * dx);
  u_x_mins_null = (u_plus_null_mins_null - u_mins_null_mins_null) ./ (2 * dx);
  u_x_null_plus = (u_plus_null_null_plus - u_mins_null_null_plus) ./ (2 * dx);
  u_x_null_mins = (u_plus_null_null_mins - u_mins_null_null_mins) ./ (2 * dx);
  
  u_y_plus_null = (u_null_plus_plus_null - u_null_mins_plus_null) ./ (2 * dy);
  u_y_mins_null = (u_null_plus_mins_null - u_null_mins_mins_null) ./ (2 * dy);
  u_y_null_plus = (u_null_plus_null_plus - u_null_mins_null_plus) ./ (2 * dy);
  u_y_null_mins = (u_null_plus_null_mins - u_null_mins_null_mins) ./ (2 * dy);
  
  # g
  g = (1 + u_x .^ 2 + u_y .^ 2) .^ -0.5;
  g_plus_null = (1 + u_x_plus_null .^ 2 + u_y_plus_null .^ 2) .^ -0.5;
  g_mins_null = (1 + u_x_mins_null .^ 2 + u_y_mins_null .^ 2) .^ -0.5;
  g_null_plus = (1 + u_x_null_plus .^ 2 + u_y_null_plus .^ 2) .^ -0.5;
  g_null_mins = (1 + u_x_null_mins .^ 2 + u_y_null_mins .^ 2) .^ -0.5;
  
  g_x = (g_plus_null - g_mins_null) ./ (2 * dx);
  g_y = (g_null_plus - g_null_mins) ./ (2 * dy);

  # fuggveny visszateresi ertek
  
  f = 0.25 * g_x .* u_x + 0.25 * g_y .* u_y + g .* (u_x_x + u_y_y);
  
endfunction