% 2008.04.25                                        2006037346 Lee Won Il
% Finding Roots Of Polynomial with Bairstow's Method
clear

n= input('Order of polynomial = ');
a= input('Coefficients (form of rowvector) = ');
r_0= input('Initial value r = ');
s_0=input('Initial value s = ');
es=input('es( Tolerance of error ) = ');
k=1; con=0; iter=0;
r=r_0; s=s_0; imax=200;

% n = Order of Polynomial
% a = Coefficients [a1, a2, ..., an+1]
% r_0,s_0 = Initial values 
% es = Tolerance of error
% k = Rowvector of roots
% con = Condition of inner loop. con=1 means it nearly converged to root
% Iter = Number of Iteration
% Imax = Maximym of Iteration

if n<3        % Input Argument Check
    error('Low order polynomial. n must be larger than 3');
end
if length(a)~=n+1
    error('Miswriting of a, Check it')
end

while n>2     % Outer Loop operate when Order > 2 

    while con==0    % Inner Loop operate when ear&eas are smaller than es
        
        iter = iter+1;

        b(1) = a(1);
        b(2) = a(2)+r*b(1);

        for i=3:n+1
            b(i) = a(i) + r*b(i-1) + s*b(i-2);  % Define b vector
        end
        
        c(1) = b(1);
        c(2) = b(2)+r*c(1);
        for i=3:n
            c(i) = b(i)+r*c(i-1)+s*c(i-2);  % Define c vector
        end

            d_0  = det([c(n-1) c(n-2);c(n) c(n-1)]);
            dr   = det([-b(n) c(n-2); -b(n+1) c(n-1)])/d_0; % delta r
            ds   = det([c(n-1) -b(n); c(n) -b(n+1)])/d_0; % delta s
            ear  = abs(dr/r)*100; % ea of r
            eas  = abs(ds/s)*100; % ea of s

        if  (ear<=es) & (eas<=es) % error test
            con = 1;
        else
            r    = r + dr; % next iteration r value
            s    = s + ds; % next iteration s value
        end
       
        if imax<iter
            break
        end 
    end   % Loop out when con=1

      % Using Quadratic Formula 
        x(k)=(r+sqrt(r^2+4*s))/2;  % Save root at rowvector x(k)
        k=k+1; % row + 1
        x(k)=(r-sqrt(r^2+4*s))/2;  % Save another root at x(k+1)
        k=k+1;

        n = n-2; 
        
        if n>2
           clear a
           a = b;
           r=r_0; s=s_0;  % reset to orgin value that we selected.
           con = 0; % reset con value
        end
    
     

end     % Loop Out when n<2;

switch n
    case 1
        x(k)=-b(2)/b(1);
    case 2
    x(k)=(-b(2)+sqrt(b(2)^2+4*b(1)*b(3)))/(2*b(1));
    k=k+1;
    x(k)=(-b(2)-sqrt(b(2)^2+4*b(1)*b(3)))/(2*b(1));
    k=k+1;
end

disp('   [Roots of Polynomial] ')
disp(x') % Transpose x : Display form of Column vector