function defl_pol_newt
% Solve P(x)=0 by method of deflation
% where P(x) is a polynomial a_{n} x^{n}+a_{n-1}x^{n-1}+ ... + a_{1}x^{1}+a_{0}
% a - is the array of the coefficients: since matlab does not allow array index to be zero, we have to use  
% a_{0}=a(1), a_{1}=a(2), ... , a_{n} = a(n+1), in particular, 
%a=[-10, 2, 4, -7, 3];        % setting up a polynomial
% corresponds to a polynomial P(x)=3x^4-7x^3+4x^2+2x-10
%a=[-136 -85 -9 5 1];
%%
a=[ 1 0 1]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = size(a); % this is to count how many coefficients came in
n = n(2);    % for some reason size of array [1, 2, 3, 4, 5, 6] is an array itself [1, 6] (1 row, 6 columns)
             % since we actually need 6 , we take the second element in size(a);  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
roots=[]     % no roots know in the beginning 
%%%%%%%%%%%% METHOD OF DEFLATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
for k=1:1:n-1   
  % solve for the root of the current polynomial by NEWTON's method, stop when subsequent approximations differ in .00001  % NEWTON'S
  p_old=2;                                                                                                                 % NEWTON'S
  p=1+sqrt(-1);                                                                                                                     % NEWTON'S  
  while abs(p-p_old)>0.00001                                                                                               % NEWTON'S
     p_old=p;                   % remember the old approximation                                                           % NEWTON'S
     [b0, b]=horners_method(a,p);     % use Horner's method to calculate P(p)                                              % NEWTON'S   
     [c0, c]=horners_method(b,p);     % use Horner's method to calculate P'(p)=Q(p)                                        % NEWTON'S  
     p=p-b0/c0;             % Newton's method: p=p-f(p)/f'(p)                                                              % NEWTON'S 
  end                                                                                                                      % NEWTON'S
roots = [roots, p];             %add the found approximate root p to the array of roots ,                             
[a0, a] = horners_method(a,p);  %calculate the coefficients of the `defalted' polynomial,  
end 
roots


% This function implemets Horner's method
% a -- (input variable) array of coefficients of the polynomial P(x)
% x0 -- (input variable) point at which Horner's matho must be implemented 
% 
% b -- (otuput variable) array of coefficienct of the polinomial Q(x): 
%       P(x)=(x-x0)Q(x)+b0,  b0 -- is the first element in the array b
function [b0, b] = horners_method(a, x0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = size(a); % this is to count how many coefficients came in
n = n(2);    % for some reason size of array [1, 2, 3, 4, 5, 6] is an array itself [1, 6] (1 row, 6 columns)
             % since we actually need 6 , we take the second element in size(a);  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b=[a(n)];                     % initialize b as array of just one element b_n=a_n  
for k=n-1:-1:2
    b=[a(k) + b(1)*x0, b] ; % ***see the explanations below
end     
if n>1        % chech for the trivial case when  the polynomial is just a constant (one coefficient) 
    b0 = a(1) + b(1)*x0;    % last coefficient
else 
    b0=a(n); b=[];          % there is only one last coefficient, polinomial Q(x) does not exist 
end     
    %***
    %a(i)+b(n-1)*x0 calculates the new element according to the formula b_{n-1}=a_{n-1}+b_{n}x0
    % (see below that b(1) is the right thing for b_{n}) 
    % [ a(i)+b(n-1)*x0, b] -- adds the new number to array b, notice, it adds new stuff the the array from the left
    % (if b=[3,4,5,6], then [2, b] works like this = [2, [3, 4, 5, 6]]= [2, 3, 4, 5, 6])
    % so, the array b is being populated from left, the newly determined coefficient is always the first element of b: b(1)