function[unew,ushort]=Firstiterate...
    (n,nbig,ce, ht, un, unlast,xi1,xi2, a,mu,is, we,wn,cc,pp,v,le)
%   given the solution un at time t_(i) and preceeding times - array un,
% computes the solution unew at t_(i+1) -vector unew,
%   using the explicit Euler method . The integral
%   is obtained by a gaussian quadrature- function integ.
% un(ti+1)= un(ti) + ht ( I(ti,x) - un(i) + Integrate( Ker(x). S(u(x,ti - tau))) ), where
% ht- discretization step in time; Ker - interaction kernel; 
% S , I - given functions 
% The vector is first computed at a grid of dimension nXn
% an then projected into a grid of dimension nbigXnbig.
% v is the propagation speed
% unew is the new layer of the solution and ushort
% is the interpolation of unew at the Chebyshev points , dimension n^2
 %%%% is is number of the last iteration in time
 ti=(is-1)*ht;
 % un1 is the last row of un, with dimension nbig ^2 X nbig^2
 for i=1:nbig^2
     un1(i) =un(is,i);
 end
%%%% unlast is the last row of un, interpolated at the Chebyshev nodes
 %%%% dimension n ^2 X n^2
 for i=1:n
     si=cos((2*i-1)*pi/(2*n));
    for j=1:n
            sj=cos((2*j-1)*pi/(2*n));
    b1(i,j)=integ(nbig,si,sj,un,un1,xi1,xi2,a,mu,we,wn,v,is,ht)-...
    unlast((i-1)*n+j);
    end
 end
%%%%% obtain ushort from b1
for i=1:n
for j=1:n
    ind=(i-1)*n+j;
    ushort(ind)=ht/ce*b1(i,j)+unlast(ind);
end
end
%%%%%%%%%%%%%%%%%  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;        
%%%%%%%%%%%%%%%%%%%%%  obtain new solution unew
%%%% convert to vector
for i=1:nbig
for j=1:nbig
    ind=(i-1)*nbig+j;
    bb(i,j)=ht/ce*bb(i,j)+un(is, ind);
    unew(ind)=bb(i,j);
end
end

