%This computes the values in Table 11


close('all')
clear all
i=sqrt(-1);

% Set up basic data

% number of points per unit length of [0,a]
hx=10000

% length of interval
a=1

% total number of points
Nx=hx*a+1

% specify the x variable
x=zeros(Nx,1);
for r=1:Nx
   x(r,1)=(r-1)/hx;
end;


% specify the initial function
f=zeros(Nx,1);
for r=1:Nx
   f(r,1)=x(r,1);
     % change end values for trapezium rule integration
   f(1,1)=f(1,1)/sqrt(2);
   f(Nx,1)=f(Nx,1)/sqrt(2);
end;

%normalize the function
f=f/norm(f);
  
    
%length of the interval S starting at 1 is called Ns
M=10
Ns=2*M+1  


% specify the un-normalized pseudo-eigenvectors
s=(1:Ns)-ones(1,Ns)*(Ns+1)/2;
%   u is the overcompleteness factor
Nu=20
values=zeros(Nu,4);
for nu=1:Nu
    nu;
    u=(0.7+nu/10);
    e=zeros(Nx,Ns);
    for rx=1:Nx          
       e(rx,:)=exp(2*pi*i*s*x(rx,1)/u  );
       e(1,:)=e(1,:)/sqrt(2);
       e(Nx,:)=e(Nx,:)/sqrt(2);
    end;

% normalize the columns of  e to unit length
for rs=1:Ns
   enorm=norm(e(:,rs));
   e(:,rs)=enorm^(-1)*e(:,rs);
end;

%compute the sequence phi for a weight delta
delta=10^(-10);
ftilde=[f;zeros(Ns,1)];
etilde=[e;delta*eye(Ns,Ns)];

phi=etilde\ftilde;


b=norm(f-e*phi);
c1=norm(phi,1);
c2=norm(phi);

values(nu,:)=[u b c1 c2];

end;

format short g

values

hold on;
plot(1:Ns,abs(phi));