%Second-order Chebyshev iteration method
%No preconditioning.
%Parameters: rhs-right-hand side; x-the solution; r-residual (Ax-b);
%            x0-initial guess; n-size of the problem
% a,b - Chebychev method parameters:
%     --- a<=min eigvalue; b>=max eigvalue;
%     --- a<=left (smaller) focus; b>= right (bigger) focus;
%  it0  - iterations performed; 
function [x]=Chebyshev(A,rhs,x0,a,b,it0)
[n,col]=size(A);
%---------------------- initialization
x(:,1)=x0;
r(:,1)=rhs-A*x(:,1);
apb=(a+b)/2;
amb=((b-a)/4)^2;
beta(1)=4/(a+b);
x(:,2)=x(:,1)+0.5*beta(1)*r(:,1);
r(:,2)=rhs-A*x(:,2);
qq=apb-amb*beta(1);
beta(2)=1/qq;
alfa(2)=apb*beta(2);
it=2;
    while (it < it0+1),
        x(:,it+1)=alfa(it)*x(:,it)+(1-alfa(it))*x(:,it-1)+beta(it)*r(:,it);
        it=it+1;
        r(:,it)=rhs-A*x(:,it);
        qq=apb-amb*beta(it-1);
        beta(it)=1/qq;
        alfa(it)=apb*beta(it);
%        pause
    end

   
    % Limit value of beta: beta_lim=4/((sqrt(a)+sqrt(b))^2);