function grad_f = diff_fin_centre(f, x, epsilon)
% calcule le gradient de 'f' au(x) point(s) 'x' par
% différences finies centrées, en utilisant un pas 'epsilon'

% x est une matrice n x m,
% n : dimension
% m : nombre de points où l'on veut calculer le gradient,
%     un point par colonne
[n, m] = size(x);

% le gradient de f au point x_j a la même dimension que x_j,
% et on calcule le gradient en m points.
grad_f = zeros(n,m);

for i = 1:n

    % construction de m exemplaires du vecteur i de la base canonique
    ei = zeros(n, m);
    ei(i,:) = 1;

    % calcul de la dérivée partielle selon la ième coordonnée,
    % pour tous les points simultanément.
    grad_f(i,:) = ( f(x+epsilon*ei) - f(x-epsilon*ei) ) / (2*epsilon);
end
