calculate a n dimensional jacobian matrix without the symbolic toolbox in matlab

main.m
function J = jacobian(f, x)
% Numerically compute the Jacobian matrix of function f(x)
% using the central difference method with a small step size h

% Initialize Jacobian matrix
n = numel(f(x));
m = numel(x);
J = zeros(n, m);

% Set step size for numerical differentiation
h = 1e-6;

% Compute Jacobian matrix using central difference method
for i = 1:m
    % Compute partial derivative of f with respect to xi using
    % central difference method
    xi = x(i);
    xplus = x;
    xminus = x;
    xplus(i) = xi + h/2;
    xminus(i) = xi - h/2;
    df = (f(xplus) - f(xminus))/h;
    
    % Store partial derivative in ith column of Jacobian matrix
    J(:, i) = df(:);
end
end
681 chars
28 lines

main.m
% Define function f(x) = [x1*x2; x1^2; x1 + x2; x2^3]
f = @(x) [x(1)*x(2); x(1)^2; x(1) + x(2); x(2)^3];

% Choose a value of x where to evaluate the Jacobian matrix
x = [1; 2];

% Compute the Jacobian matrix of f at x using numerical differentiation
J = jacobian(f, x)
270 chars
9 lines

main.m
J =

    2.0000    1.0000
    2.0000         0
    1.0000    1.0000
         0   12.0000
89 chars
7 lines