function [un,ushort]=Alliterates...
    (n,nbig,un,unlast,unlast1, xi1,xi2, a, mu,is,ce,ht,we,wn,cc,pp,eps,v)
% given the solution at time t(i-1) and t(i-2), and previous times
% computes the solution unew (time ti).
%  To  solve the nonlinear uses the fixed point method.
% The initial approximation un0 is obtained by the explicit Euler method.
% is is the index of the last computed iterate (at time ti-1)
% ushort is un0 at the Chebyshev nodes
% unlast and unlast1 are the two last rows of the solution at
% Chebyshev points
[un0,ushort]=...
Firstiterate(n,nbig,ce,ht,un,unlast,xi1,xi2,a,mu,is,we,wn,cc,pp,v);
     for k=1:nbig^2
         err(k)=abs(un0(k)-is*ht);
     end
   %  normerr=max(err)
it=0;
dif=1;
ti=is*ht;
%%%%%%%%%%%%%%%%  un1 is the vector with the last row of the solution
for i=1:nbig^2
     un1(i) =un(is,i);
end
un1short=ushort;
%%%% start the fixed point method
while(it<10 & dif>eps)
    it=it+1;
for i=1:n
     si=cos((2*i-1)*pi/(2*n));
    for j=1:n
         sj=cos((2*j-1)*pi/(2*n));
        ind=(i-1)*n+j;
 b1(i,j)=I(si,sj) -un1short(ind)+...
     +integ(nbig,si,sj,un,un1,xi1,xi2,a,mu,we,wn,v,is+1,ht);
    end
end
%%%%%%%%%%%%  build ushort from b1
for i=1:n
for j=1:n
    ind=(i-1)*n+j;
 ushort(ind)=(2*ht)/(3*ce)*b1(i,j)+(4*unlast(ind)-unlast1(ind))/(3*ce) ;
end
end
un1short=ushort;
%%%%%%%%%%%%%%%%%  build matrix lambda
%%%%%%%%%%%%%  cc is a matrix n X n
l=b1*cc';
la=cc*l;
%sla=size(la)
%%%%%%%%%%%%%%%%  transform matrix (polynomial interpolation)
bb=pp'*la;
%sbb=size(bb)
bb=bb*pp;  
%%%%%%%%%%%%%%%%%%%%%  add to the previous iteration
%%%% convert to vector
for i=1:nbig
for j=1:nbig
    ind=(i-1)*nbig+j;
    bb(i,j)=2*ht/(3*ce)*bb(i,j)+(4*un(is,ind)-un(is-1,ind))/(3*ce);
    unew(ind)=bb(i,j);
end
end
%%%%%%%%%%%%%%%%%% diference between iterates
dif=norm(unew-un1);
un1=unew;
%%%%%%%%%%% end of cycle while    
end
for i=1: nbig^2
    un(is+1,i)=unew(i);
end

