function dy = mfs_derive(y, dx)

# usage: dy = mfs_derive(y, dx)
#
# Input  y(:, :)   Array with function values
#                  (Rows correspond to x-values)
#        dx        x-increment
#
# Output dy(:, :)  Array with derivatives of y
#
# The function computes the first derivatives using central
# differences.
#
# -------------------------------------------------------------------

  n = columns(y);

  dy(:, n) = y(:, n) - y(:, n - 1);
  dy(:, 1) = y(:, 2) - y(:, 1);

  dy(:, 2 : n - 1) = 0.5 * (y(:, 3 : n) - y(:, 1 : n - 2));

  dy = dy / dx;

end
