function J = numerical_jacobi(f,x,eps)
xe = zeros(size(x));
xe(1) = eps;
dy = (f(x+xe)-f(x-xe))./eps./2;
J = zeros(length(dy),length(x));
J(:,1) = dy;
for i=2:length(x)
    xe(i-1) = 0;
    xe(i) = eps;
    dy = (f(x+xe)-f(x-xe))./eps./2;
    J(:,i) = dy;
end
