function  DD=ddd_tabl(FUN,dFUN,xn);
%%  function  DD=ddd_tabl(FUN,dFUN,xn);
%%
%%  THE DIVIDED DIFFERENCE TABLE for the function and 
%%  derivative values at points with multiplicity 2
%%  (each point taken twice)
%%  INPUT:
%%   xn   - the points at which the divided difference is 
%%          taken as a row vector ,   e.g. [0,1,2,3,4]
%%    FUN  - the function value data for the difference which 
%%           may be entered EITHER as a string, enclosed in ' ',
%%           with the name the m-file of the function,  OR as a  
%%           row vector of data in the same order as xn
%%   dFUN  - the derivative value data for the difference which
%%           may be entered EITHER as a string, enclosed in ' ',
%%           with the name the m-file of the derivative function, 
%%           OR as a row vector of data in the same order as xn
%%  OUTPUT:
%%    DD  - matrix containing the divided difference table, 
%%          with xn taken twice as first column, FUN at xn 
%%          xn taken twice as second column, and dFUN at xn 
%%          occupying the even rows of column 3 ... the k-th 
%%          order divided difference is in the (k-2)nd column

n=length(xn)-1;  % determine the number of points
yn(1:2:2*n+1)=xn;  % double up the points
yn(2:2:2*n+2)=xn;

DD(:,1)=yn';   % points as first column in the table, each repeated

if isstr(FUN)==1; % test whether string or data
   DD(:,2)=feval(FUN,yn');  % 2nd column are the function values
else
   DD(1:2:2*n+1,2)=FUN.'; % 2nd column as given data
   DD(2:2:2*n+2,2)=FUN';
end;

for i=2:2*n+2;       %% an adjustment is made for indexing in MATLAB
  for j=3:i+1;
     if (i/2==floor(i/2) & j==3)
	if isstr(dFUN)==1;
           DD(i,j)=feval(dFUN,xn(i/2));
	   % derivative values as 3nd column at the
	   %  second occurence of each point
        else
           DD(i,j)=dFUN(i/2);  
	end;
     else
        DD(i,j)=(DD(i,j-1)-DD(i-1,j-1))/(yn(i)-yn(i-j+2));
     end;
  end;
end;
