function [A,B,C,D] = central_difference(x0,u0,par)
format long
dx_value = 1e-7;
du_value = 1e-7;

Ns=3; %Number of states
Nm=8; %Number of measurements
Nu=6; %Number of inputs

A = zeros(Ns,Ns);
B = zeros(Ns,Nu);
C = zeros(Nm,Ns);
D = zeros(Nm,Nu);

dx= zeros(Ns,1);
for i =1:Ns
    dx(i)=dx_value;
    x1 = x0 + dx;
    x2= x0 - dx;
    xdot1=v1_gaslift_model([],x1,u0,'derivatives',par);
    xdot2=v1_gaslift_model([],x2,u0,'derivatives',par);
    
    y1 = v1_gaslift_model([],x1,u0,'measurements',par);
    y2 = v1_gaslift_model([],x2,u0,'measurements',par);

    A(:,i) = (xdot1-xdot2)/(2*dx(i));
    C(:,i) = (y1-y2)/(2*dx(i));
    dx(i)=0;
end

du=zeros(Nu,1);
for i =1:Nu
    du(i)=du_value;
    u1 = u0 + du;
    u2= u0 - du;
    
    xdot1=v1_gaslift_model([],x0,u1,'derivatives',par);
    xdot2=v1_gaslift_model([],x0,u2,'derivatives',par);
   
    y1 = v1_gaslift_model([],x0,u1,'measurements',par);
    y2 = v1_gaslift_model([],x0,u2,'measurements',par);
    
    B(:,i) = (xdot1-xdot2)/(2*du(i));
    D(:,i) = (y1-y2)/(2*du(i));
    du(i)=0;
end
%Remove last column which is just transfer function from P0
%B(:,end)=[];
%D(:,end)=[];

end

