function[pochodna] = stirling(x,y,h)
    n = size(x,2);
    T=zeros([size(x,2) size(x,2)+1]);
    T(:,1) = x;
    T(:,2) = y;
    fx = 0;
    for j=3:size(x,2)+1
        for i=size(x,2)-(j-2):-1:1
            T(i,j)=T(i+1,j-1)-T(i,j-1);
        end
    end
    a=[1,-1/6,1/30];
        for i=2:2:n
            fx = fx + a(i/2)*(0.5*(T((((n-i)+1)/2),i+1)+T((((n-i)+1)/2+1),i+1)));
        end
        pochodna = fx * 1/h;
end