function dp = dempol(n,xmin,xmax,alp,verbose)
%
% n = degree of polynomial
%
% Polynomial interpolation using Neville's algorithm
%
% Procedure POLINT from Numerical Recipes
%
% Input: n equidistant points between xmin and xmax
%
% Y given by function exp(x/alp) 
%
% For various X, the process of calculation of the interpolation 
% polynomial is demonstrated
    
xval=zeros(1,n);
yval=zeros(1,n);
    
dx = (xmax-xmin)/n;
xval = (xmin:dx:xmax);
% evaluate the function for each value of vector xval
yval = funval( xval,alp );
for j=1:n+1
   fprintf('%5d:  %12.7e,  %12.7e\n',j,xval(j),yval(j));
end
x=0;
while (x > -1.e+30)
    x = input('Enter x   (end by x <= -1.e+30): ');
    if (x > -1.e+30)
        [y, dy] = polint(xval,yval,n+1,x,verbose);
        fprintf(' For x = %11.4e is y = %11.4e with estim. error= %11.4e \n',...
        x,y,abs(dy));
        fprintf(' function value  %11.4e \n',funval(x,alp));
    end
end
    
%%%%%%%%%%%%

function fv = funval(x,alp);
%
% The function to be approximated. 
%    
% Works for scalar as well as vector arguments
%
    fv = exp(x/alp);
    return;
    
%%%%%%%%%%%%

% $$$ THE FOLLOWING WAS MOVED TO SEPARATE FILE

% $$$ function [y,dy] = polint(xa,ya,n,x);
% $$$ %
% $$$ % Calculates value Y of the interpolation polynomial of order N-1
% $$$ % at point X
% $$$ %    
% $$$ % dY is the error estimate of Y
% $$$ %
% $$$ % Given function values are xa(n), ya(n)
% $$$ %
% $$$     
% $$$ c = zeros(n);
% $$$ d = zeros(n);
% $$$     
% $$$ ns  = 1;
% $$$ dif = abs(x-xa(1));
% $$$ fprintf('Step  1 \n');
% $$$ 
% $$$ for i = 1:n
% $$$     dift = abs(x-xa(i));
% $$$     if (dift<dif)
% $$$         ns  = i;
% $$$         dif = dift;
% $$$     end;
% $$$     c(i) = ya(i);
% $$$     d(i) = ya(i);
% $$$     fprintf(' c(%2d) = %12.7e,  d(%2d) = %12.7e \n',i,c(i),i,d(i));
% $$$ end;
% $$$ y = ya(ns);
% $$$ fprintf(' ns = %2d,  value y in 1. step is %12.7e\n',ns,y);
% $$$ ns = ns-1;
% $$$ for m=1:n-1
% $$$     pause;
% $$$     fprintf('Step %2d \n')
% $$$     for i=1:n-m
% $$$         ho  = xa(i)  -x;
% $$$         hp  = xa(i+m)-x;
% $$$         w   = c(i+1)-d(i);
% $$$         den = ho-hp;
% $$$         if (den==0.0)
% $$$             fprintf('pause in routine POLINT');
% $$$             pause;
% $$$         end;
% $$$         den  = w/den;
% $$$         d(i) = hp*den;
% $$$         c(i) = ho*den;
% $$$         fprintf(' c(%2d) = %12.7e,  d(%2d) = %12.7e \n',i,c(i),i,d(i));
% $$$     end;
% $$$ 
% $$$     if (2*ns < n-m) 
% $$$         dy = c(ns+1);
% $$$         fprintf(' dy = c(%2d) = %12.7e \n',ns+1,dy);
% $$$     else
% $$$         dy = d(ns);
% $$$         fprintf(' dy = d(%2d) = %12.7e \n',ns,dy);
% $$$         ns = ns-1;
% $$$     end;
% $$$     y = y+dy;
% $$$     fprintf(' Value y in step %2d is %12.7e \n',m+1,y);
% $$$ end;
% $$$ return;
    
