gauss forward interpolation method

I wrote the following code for gauss interpolation

function [ yval ] = gauss_p( xd,yd,xp )
n=length(xd);
if(length(yd)==n) 
tbl=yd';
for j=2:n
    for i=1:n-j+1
        tbl(i,j)=tbl(i+1,j-1)-tbl(i,j-1);
    end
end
tbl
h=xd(2)-xd(1);
if rem(n,2)==0
    k=n/2+1;
else
    k=n/2+0.5;
end
p=(xp-xd(k))/h
pt=cumprod([1,p-(0:n-3)])
dt=[tbl(k,1) tbl(k,2) tbl(k-1,3)+tbl(k-1,4) tbl(k-1,4)+tbl(k-1,5) tbl(k-1,5)+tbl(k-1,6) tbl(k-1,6)+tbl(k-1,7)]./factorial(0:n-2)
yval=sum(pt.*dt)
end

Suggest me to simplify the calculation of 'dt' vector subject to 'n' values

Answer 1

[Torsten]

dt = [tbl(k,1), tbl(k,2), tbl(k-1,3:n-1)-tbl(k-1,4:n)]./factorial(0:n-2)

?

Best wishes

Torsten.