function [ un] = mainfunction(n,nbig,nt, ht,eps,mu,v)
%this is the main program
%   it builds the solution at any moment nt.
% INPUT: n -number of Chebyshev nodes (small dimension); nbig - number of
% Gaussian nodes (big dimension); nt- number of steps in time; 
%eps - tolerance for  iterative process at each step; mu -parameter for
% sigmoidal function; v - propagation velocity
%  parameters for the connectivity kernel:
xi1=0.5; 
xi2=1;
a=0;
% pop. dynamics
ce=1; 
% 1) Computing weights and nodes
we=weights(nbig);
wn=wnodes(nbig);
% 2) compute matrix P  
  pp=pmatrix(n,nbig);
   %3) matrix for cos. transform
   cc=comatrix(n);
 %4) Initialize array for solution
 % nt is the number of steps in time
  un=zeros(nt+1,nbig^2);
  % 5) set initial approximation
 for i=1:nbig
     for j=1:nbig
         ind=(i-1)*nbig+j;
     un(1,ind)=1e-2;
     end
 end
        for i=1:n^2
            for j=1:n^2
         ind=(i-1)*nbig+j;
              un0(ind)=1e-2;
            end
        end
%6) Compute solution at instant t1=ht
[unew,ushort]=...
Firstiterate(n,nbig,ce, ht,un,un0,xi1,xi2,a,mu,1, we,wn,cc,pp,v);
%%%%%%%%%%%%%%% ushort is the solution at the Chebyshev nodes
 %%%%%%%%%%%%% 7) store unew in solution array   
for i=1:nbig^2
      un(2,i)=unew(i);
end
  %%%%%%%%% 8) initialize multistep method
  unlast=ushort;
  unlast1=un0;
 %%%%%%%%%%%%%%%%%%  9) Computing solution from time t2 
 for i=2:nt
     % i is the index of the solution at time t_(i-1)
[un,ushort]=Alliterates...
(n,nbig,un,unlast,unlast1,xi1,xi2,a,mu,i,ce,ht,we,wn,cc,pp,eps,v);
%%%%%% update arrays
     unlast1=unlast;
     unlast=ushort;
' t='
    i*ht
     for k=1:nbig^2
      ulong(k)=un(i+1,k);
  end
   umax=max(ulong)
   umin=min(ulong)
 end
end