function T = difference_table(Y)
n = length(Y);
%Construct empty forward difference table
T = zeros(n, n);
%Fill first column
T(:,1) = Y;
%Fill remaining columns
for j = 2:n % j is the column index
    for i = j:n % i is the row index
        T(i,j) = T(i,j-1) - T(i-1,j-1);
    end
end