function [A,B,C,D] = v4_new_4d_linmod_num(x0,u0,par)


A = zeros(4,4);
B = zeros(4,3);
C = zeros(4,4);
D = zeros(4,3);

dx = [0;0;0;0];
for i =1:4
    dx(i) = 1e-2;
    x = x0 + dx;
    xdot0 = v4_new_4d_model(0,x0,u0,'derivatives',par);
    xdot = v4_new_4d_model(0,x,u0,'derivatives',par);
    y0 = v4_new_4d_model(0,x0,u0,'measurements',par);
    y = v4_new_4d_model(0,x,u0,'measurements',par);
    A(:,i) = (xdot-xdot0)/dx(i);
    C(:,i) = (y(1:4)-y0(1:4))/dx(i);
    dx(i) = 0;
end

du = [0;0;0];
for i =1:3
    du(i) = 1e-2;
    u = u0 + du;
    xdot0 = v4_new_4d_model(0,x0,u0,'derivatives',par);
    xdot = v4_new_4d_model(0,x0,u,'derivatives',par);
    y0 = v4_new_4d_model(0,x0,u0,'measurements',par);
    y = v4_new_4d_model(0,x0,u,'measurements',par);
    B(:,i) = (xdot-xdot0)/du(i);
    D(:,i) = (y(1:4)-y0(1:4))/du(i);
    du(i) = 0;
end

end