function out=Jaco(x)
% jacobian calls f.m
% f.m gives column vectors of the different equations to calculate.
% inputs: vector of x values 
% outputs: out=jacobian using centered difference

h=0.0001; % sets the interval on which we want to calculate centered difference.
n=length(x); % gives the dimension of the jacobian

for j=1:n
    e=[zeros(1,j-1),h,zeros(1,n-j)]; % makes a vector [0,0,...,0,h,0,0,0,...,0] h is at the j_th position
    xjp=x+e; % gives a vector similar to e, but has the x values instead of zeros ie. [x1,x2,...,h,...xn]
    xjm=x-e; % [x1,x2,...,-h,...xn]
    fjp=f(xjp); % gives the function value at f(x+h)
    fjm=f(xjm); % gives the function value at f(x-h)
    out(:,j)=(fjp-fjm)/2/h; % centered difference formula
end