% Section 3.1, problem 6(a)

format long
x=[ 0    0.25    0.5     0.75];
y=[ 1 1.64872 2.71828 4.48169];
xstar = 0.43;

n = length(x)-1; 
Q = zeros(n+1);
Q(:,1) = y';
R = Q;

for i = 1:n
  for j = 1:i
    Q(i+1,j+1) = ( (xstar-x(i-j+1))*Q(i+1,j) - (xstar-x(i+1))*Q(i,j) ) ...
                /( x(i+1)-x(i-j+1) );
  end
end        
Q(n+1,n+1)

for j = 1:n
  for i = n:-1:j
    R(i+1,j+1) = ( (xstar-x(i-j+1))*R(i+1,j) - (xstar-x(i+1))*R(i,j) ) ...
                /( x(i+1)-x(i-j+1) );
  end
end        
R(n+1,n+1)