function A8

%number of interp. nodes -1
n=7;
%prepare plot
figure;
%do plain equidistant polynomial interpolation
% interpolation points
I=[1:(n+1)]-1;
p=I/n;
% function to be interpolated
y=f(p);
% number of evaluation points
N=100;
% evaluation points
x=linspace(0,1,N);
%prepare polynomial
s=zeros(1,length(x)); 
%combine Lagrange polynomials
for i=1:length(p) s=s+y(i)*lag(i,p,x); end; 

%do ENO type reconstruction
%order
o=4;
%add m-1 more points for boundary effects
m=o+2;
p=[[-m+1:-1]*p(2),p,1+[1:m-1]*p(2)];
%calculate function values
y=f(p);
%prepare result vectors
X=[]; Y=[];
Y2=[];
for i=m:m+n
%take 10 nodes per intervall for plotting
xs=linspace((p(i)+p(i-1))/2,(p(i)+p(i+1))/2,10);
%select the interpolation nodes for the interval
[nodes,coeff]=select(y,p,i,o);
%evaluate Newton form
P=newton_ev(p(nodes),coeff,xs);
%evaluate Newton form for fixed interpolation nodes (left, middle, right)
Q=newton_ev([p(i),p(i-1),p(i+1)],[y(i),divdif([y(i),y(i-1)],[p(i),p(i-1)]),divdif([y(i),y(i-1),y(i+1)],[p(i),p(i-1),p(i+1)])],xs);
%combine with data for previous intervals
X=[X,xs];
Y=[Y,P];
Y2=[Y2,Q];
end;
%now plot everything
plot(p,y,'bo',X,Y,x,f(x),'g',X,Y2,'c',x,s,'y');
axis([0,1,-4,4]);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [y]=f(x)
%function to be interpolated

I=find(x<0.5);
y(I)=exp(3*x(I))-1;
I=find(x>=0.5);
y(I)=sin(6*x(I)-3)-3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [nodes,coeff]=select(y,x,i,o)

%number of interpolation points (order+1)
n=o+1;
%indices of the selected points
nodes=zeros(1,n);
%coefficients of corresponding Newton form
coeff=zeros(1,n);
nodes(1)=i;
%prepare tableau
T=zeros(n,n);
T(1,1)=y(i);
coeff(1)=T(1,1);
%fill tableau line by line
for i=2:n
%prepare two tableaus
H1=T;
H2=T;
%find first and last node
first=min(nodes(1:i-1))-1;
last=max(nodes(1:i-1))+1;
%put corresponding function values in last line of the tableau
H1(i,1)=y(first);
H2(i,1)=y(last);
%coordinates corresponding to selected nodes
X=x(nodes(1:i-1));
%calculate last lines of tableaus
for j=2:i
H1(i,j)=(H1(i,j-1)-H1(i-1,j-1))/(x(first)-X(i-j+1));
H2(i,j)=(H2(i,j-1)-H2(i-1,j-1))/(x(last)-X(i-j+1));
end;
a=abs(H1(i,i));
b=abs(H2(i,i));
%throw the bad tableau and keep the good one
if a<=b T=H1; nodes(i)=first; else T=H2; nodes(i)=last; end;
%remember the divided difference
coeff(i)=T(i,i);
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [d]=divdif(y,x)

%number of interpolation points
n=length(x);
%prepare tableau
T=zeros(n,n);
%fill tableau line by line
for i=1:n
T(i,1)=y(i);
for j=2:i
T(i,j)=(T(i,j-1)-T(i-1,j-1))/(x(i)-x(i-j+1));
end;
end;
d=T(n,n);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [y]=lag(i,p,x)
% evaluate Lagrange polynomial i defined by 
% l(i,p,x(j))=delta(i,j)
% at the position x and return value
% y=l(i,p,x)

%number of points
m=length(p);
%get vector of differences
q=[p(1:i-1),p(i+1:m)];
q=p(i)-q;
% and invert it
q=1./q;
% multiply to get coefficient 
c2=prod(q);
% form of vector x
[ix,jx]=size(x);
nx=ix*jx;
%make it a long row vector and replicate m-1 times
X=repmat(x(:)',m-1,1);
% get matrix of differences to x
q=repmat([p(1:i-1),p(i+1:m)]',1,nx);
q=X-q;
% calculate Lagrange polynomial
if m>2 psi=(c2*prod(q)); else psi=c2*q; end;
% bring back to original form
y=reshape(psi(:),ix,jx);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [y]=newton_ev(p,c,x);
%evaluate the Newton form with point coordinates p and 
%coefficients c at the points x

%number of points
m=length(p);
% form of vector x
[ix,jx]=size(x);
nx=ix*jx;
%make it a long row vector and replicate m-1 times
X=repmat(x(:)',m-1,1);
% get matrix of differences to x
q=repmat(p(1:m-1)',1,nx);
q=X-q;
%values of Newton basis polynomials
for i=2:m-1 
q(i,:)=q(i,:).*q(i-1,:);
end;
%multiply with coefficients
y=c(1)+c(2:end)*q;
