How to calculate Jacobian matrix through numerical analysis.

Hello. I would like to obtain the Jacobian matrix of the aircraft and implement the extended Kalman filter. I want to estimate the state variable after n seconds through the Extended Kalman Filter. I implemented the following formula.

xdot is , d_xdot is , delta is.

function [J,fdx,fx,dx,deltafx] = Jacobian(xdot,d_xdot,delta)

fdx=zeros(9,1);
fx=zeros(9,1);
dx=zeros(9,1);
J=zeros(9);
deltafx=zeros(9);
fdx(1) = d_xdot(1);
fdx(2) = d_xdot(14);
fdx(3) = d_xdot(27);
fdx(4) = d_xdot(40);
fdx(5) = d_xdot(53);
fdx(6) = d_xdot(66);
fdx(7) = d_xdot(82);
fdx(8) = d_xdot(95);
fdx(9) = d_xdot(108);

fx(1)=xdot(1);
fx(2)=xdot(2);
fx(3)=xdot(3);
fx(4)=xdot(4);
fx(5)=xdot(5);
fx(6)=xdot(6);
fx(7)=xdot(10);
fx(8)=xdot(11);
fx(9)=xdot(12);

dx(1)=delta(1);
dx(2)=delta(2);
dx(3)=delta(3);
dx(4)=delta(4);
dx(5)=delta(5);
dx(6)=delta(6);
dx(7)=delta(10);
dx(8)=delta(11);
dx(9)=delta(12);

for i=1:9
    for j=1:9
        J(i,j)=(fdx(i)-fx(i))/dx(j);
    end
end

end

'd_xdot' is calculated as follows.

I would like to ask for advice if this is correct.