% This function calculates the Jacobian matrix of a function ('fun') at point x0 
% This prograpm mimics the procedure "gradp" in the GAUSS library

function y=gradp(fun,x0);
 

    f0 = fun(x0);
    [r1, c1] = size(f0);
    n = r1;
    [r2, c2] = size(x0);
    k = r2;
    grdd = zeros(n,k);

% Computation of stepsize (dh) for gradient 

    ax0 = abs(x0);
    if x0 ~= 0;
        dax0 = x0./ax0;
    else;
        dax0 = 1;
    end;
    
    dh1 = ones(k,2);
    dh1(1:k,1)=ax0;
    dh = (1.e-8)*max(dh1')'.*dax0;
    xdh = x0+dh;
    dh = xdh-x0;    % This increases precision slightly 
    
    
    arg=zeros(k,k);
    
    i=1;
    while i <=k;
    
        arg(:,i)=x0;
      i = i+1;  
    end;

    i=1;
    while i <=k;
    
        arg(i,i)=xdh(i);
     i = i+1;   
    end;
    
    i = 1;
    while i <=k;
        grdd(:,i) = fun(arg(:,i));
        i = i+1;
    end;

    grdd = (grdd-f0)./(dh');

    y=grdd';