interpolat ( func , xmin , xmax , ns , n , es , maxit )

function xr = interpolat ( func , xmin , xmax , ns , n , es , maxit )
% needs GaussSiedel.m
% finds brackets of x that contain sign changes of a function on 
% the interval [ xmin , xmax ] and uses interpolation method to find the
% root of a function
% xr = interpolat ( func , xmin , xmax , ns , n , es , maxit )
%
% arguments: ( input )
%  func - name of function
%
%  xmin - scalar ( double ) - defines the lower guesses of the root
%
%  xmax - scalar ( double ) - defines the upper guesses of the root
%
%  ns - scalar ( integer ) - ( optional ) number of sections minus one or number of
%  subintervals along x axis used to search for brackets
%
%  n - scalar ( integer ) - degree of polynomial plus one
%
%  es - scalar ( double ) - ( optional ) stopping criterion
%
%  maxit - scalar ( integer ) - ( optional ) maximum allowable iterations
%
% arguments: ( returned )
%   xr - real roots
%
% Example: 
%  fp=inline ( '(-12)+(-21.*x)+(8*x.^2)-(2.75*x.^3)','x');
%  A=interpolat(fp,2,-1,1000,3,0.001,100)
% A = 
%    -0.47438852171570 

% check for simple errors
% if necessary, assign default values
if nargin > 7
  error 'less than 8 arguments must be provided'
end
if nargin < 7 %if maxit blank set to 50
     maxit = 50 ;
end
if nargin < 6 %if es blank set to 0.001  
    es = 0.001 ;
end
if nargin < 5 %if n blank set to 10
    n = 10 ;
end
if nargin < 4 %if ns blank set to 50
    ns = 50;
end
if ( length ( xmin ) ~= 1 ) || ( ~isfloat ( xmin ) )
  error 'xmin must be a scalar double variable'
elseif ( length ( xmax ) ~= 1 ) || ( ~isfloat ( xmax ) )
  error 'xmax must be a scalar double variable'
elseif ( length ( ns ) ~= 1 ) || ( ~isinteger ( ns ) ) || ( ns < 1 )
  error 'ns must be a scalar positive integer variable'
elseif ( length ( n ) ~= 1 ) || ( ~isinteger ( n ) ) || ( n < 1 )
  error 'n must be a scalar positive integer variable'
elseif ( length ( es ) ~= 1 ) || ( ~isfloat ( es ) )
  error 'es must be a scalar positive double variable'
elseif ( length ( maxit ) ~= 1 ) || ( ~isinteger ( maxit ) ) || ( maxit < 1 )
  error 'maxit must be a scalar positive integer variable'
end
if nargin < 3
  error 'more than 2 arguments must be provided'
end

% Incremental search
x = linspace ( xmin , xmax , ns ) ; % Devide the interval into ns - 1 section
f = func ( x ) ; % Evaluates the value of function in the points
nb = 0 ;
xb = [ ] ; % xb is null unless sign change detected
for k = 1 : length ( x ) - 1% k is the numerator of for loop
   if sign ( f ( k ) ) ~= sign ( f ( k + 1 ) ) % Checks for sign change
    nb = nb + 1 ;
    xb ( nb , 1 ) = x ( k ) ;
    xb ( nb , 2 ) = x ( k + 1 ) ;
  end % End of if loop
end % End of for loop
[ l , w ] = size ( xb ) ;
xr = [ ] ;
iter = [ ] ;
xrold = [ ] ;
A = [ ] ;
r = [ ] ;
ea = [ ] ;
% Polynomial interpolation to find roots
for i = 1 : l % l is the numerator of for loop
    xr ( i ) = xb ( i , 1 ) ;
    iter ( i ) = 0 ;
    while ( 1 )
        xrold ( i ) = xr ( i ) ; 
        xt = linspace ( xb ( i , 1 ) , xb ( i , 2 ) , n ) ;
        g = func ( xt ) ;
        for j = 1 : n % j is the numerator of for loop
            for l = 1 : n % l is the numeratot of for loop
                A(l,j)=(xt(l))^(j-1) ; % A is the generator of polynomial matrix
            end % End of for loop
        end % End of for loop
        A = fliplr(A); % Makes increasing A matrix. 
        X = GaussSeidel(A,(g')); % Evaluates the eigenvalues of matrix A.
        roo=roots(X); % Evaluates roots of polynomial.
        nm=1; % Index of roots.
        for nu=1:length(roo) % nu is the numerator of for loop
            if imag (roo(nu))==0 % Selects the real roots.
                r(nm)=roo(nu);
                nm=nm+1;
            end % End of if loop
        end % End of for loop
        ma=r(1);
        for mx = 1:length(r) % mx is the numerator of for loop
            if abs(func(r(mx)))<abs(func(ma)) % Finds the nearest value to root.
                ma=r(mx);
            end % End of if loop
        end % End of for loop
        xr(i)=ma;
        test = func(xb(i,1))*func(xr(i)); % Tests if the root is in the interval.
        if test < 0
            xb(i,2) = xr(i);
        elseif test > 0
            xb(i,1) = xr(i);
        else
            ea(i) = 0;
        end % End of if loop
        if xr ( i ) ~= 0 % Counts relative error percentage.
            ea ( i ) = abs ( (abs(func( xr ( i ))) - abs(func(xrold ( i ))) ) /abs(func( xr ( i ))))  * 100 ; 
        end % End of if loop
        iter ( i ) = iter ( i ) + 1;
        if ea ( i ) <= es || iter ( i ) >= maxit
            break % If the conditions take place terminates while loop
        end % End of if loop
    end % End of while loop
end % End of for loop
if isempty ( xb )    % Display that no brackets were found
    disp ( ' no roots found ' )
    disp( ' check interval or increase ns ' )
end % End of if loop
% With special thanks to John D'Errico
% By Ali Mohammad Razeghi