% id=0 : y=Av
% id=1 : y=p(A)Av = v-q(A)v
% id=2 : y=p(a)v
function y=cheb_fun(v,id)    
global A Theta Delta k_cheby
if id==0
   y=A*v; return
elseif id==1
   v=A*v; 
end
sigma_1=Theta/Delta;
 rho_0= 1/sigma_1; 
  rho = 1/(2*sigma_1-rho_0);       %%-------------- Repeating one step
  x_0 = v/Theta;
      if k_cheby< 0, y=v;   return; end %% EXIT No polynomials
      if k_cheby==0, y=x_0; return; end %% EXIT Degree=0 (Constant)
    w = A*v; 
  x_k = 2*rho/Delta*( 2*v - w/Theta); 
  rho_0 = rho;

for j=2:k_cheby 
  rho = 1/(2*sigma_1-rho_0);
  w=A*x_k;
  w=rho * ( 2*(sigma_1*x_k - w/Delta) - rho_0*x_0 + 2*v/Delta);
  x_0=x_k; x_k=w;   % p_{k}v guys
  rho_0=rho;
end
  y = x_k  ;
return %=======================================================
