clear all
%eigs tis NXN proseggisis tou -d^2/dx^2 in domain [0,pi]

Nf=20;
N=Nf-2;
I=eye(N);

%Aytos einai arithmos iterations. To programma einai prvtogono den vazei
%kritiria
Nstep=100;
xf=linspace(0,pi,Nf)';
x=xf(2:N+1);

delta=x(2)-x(1);

D2=diag(2/delta^2*ones(1,N))-diag(1/delta^2*ones(1,N-1),-1)-diag(1/delta^2*ones(1,N-1),1);

[U,S]=eig(D2);

Sd=diag(S);[sd,ind]=sort(real(Sd));
Sd=Sd(ind);U=U(:,ind);
Sd(1:5)


leig=linspace(1,N,N).^2;

%ypologismos toy antistrofou
invD2=I/D2;

psi=randn(N,1);

psi=psi/norm(psi);
psin=psi;


%it einai h taxi tis idiotimis
for it=1:18;
    
psi=randn(N,1);

psi=psi/norm(psi);
if it==1;
psin=psi; 
else;
    %edo ginetai h provoli sto katheto xoro
psi=psi-psieig(:,1:it-1)*(psieig(:,1:it-1)'*psi);
psin=psi/norm(psi);
end;


for ia=1:Nstep;
    psin=invD2*psin;
    
    if it==1;
    psin=psin/norm(psin);
    else;
    psin=psin-psieig(:,1:it-1)*(psieig(:,1:it-1)'*psin);
    %psin'*psieig(:,1)
    %pause
    psin=psin/norm(psin);
    end; 
end;
    fg=invD2*psin./psin;
    lambda(it)=mean(fg);
    psieig(:,it)=psin;
    figure(89)
    plot(x,psieig)
    drawnow
    
end;
    %eigs of D2:
    lambda=1./lambda;
    
    
    figure(68)
semilogy(sqrt(leig),leig,sqrt(leig),Sd,'or',sqrt(leig),lambda,'+k')
xlabel('order nunber');ylabel('\lambda')
title('solid: n^2 , red circles: eigs of 20 by 20 matrrix, black +: eigs as obtained by iteration')