Procedure TQLI(var D:array of real; E:array of real; N:integer;var Z:matrx2);
Label 1,2;
Var
    I,L,M,K,ITER:integer;  G,R,C,F,S,P,DD,ZZ,ZZZ,B:real;
begin
    If N > 1 Then
    begin
        For I:=2 To N do
            E[I - 1]:=E[I];
        E[N]:=0; 
        For L:=1 To N do
        begin
            ITER:=0;
1:          For M:=L To N - 1 do
            begin
                DD:=Abs(D[M]) + Abs(D[M + 1]);
                If Abs(E[M]) + DD = DD Then GoTo 2;
            end; 
            M:=N;
2:          If M <> L Then
            begin
                If ITER = 30 Then showMessage(' too many iterations ');
                ITER:=ITER + 1;
                G:=(D[L + 1] - D[L]) / (2  * E[L]);
                R:=Sqrt(Sqr(G) + 1 );
                If G >= 0 then
                    ZZ:=1
                Else
                    ZZ:=-1;
                If G >= 0 then
                    ZZZ:=1
                Else
                    ZZZ:=-1;
                G:=D[M] - D[L] + E[L] / (G + ZZZ * ZZ);
                S:=1;
                C:=1;
                P:=0; 
                For I:=M - 1 DownTo L do
                begin
                    F:=S * E[I];
                    B:=C * E[I];
                    If Abs(F) >= Abs(G) Then
                    begin
                        C:=G / F;
                        R:=Sqrt(Sqr(C) + 1 );
                        E[I + 1]:=F * R;
                        S:=1 / R;
                        C:=C * S;
                    end
                    Else
                    begin
                        S:=F / G;
                        R:=Sqrt(Sqr(S) + 1 );
                        E[I + 1]:=G * R;
                        C:=1 / R;
                        S:=S * C;
                    end;
                    G:=D[I + 1] - P;
                    R:=(D[I] - G) * S + 2  * C * B;
                    P:=S * R;
                    D[I + 1]:=G + P;
                    G:=C * R - B;
                    //Omit lines from here ...
                    For K:=1 To N do 
                    begin
                        F:=Z[K, I + 1];;
                        Z[K, I + 1]:=S * Z[K, I] + C * F;
                        Z[K, I]:=C * Z[K, I] - S * F;
                    end; 
                   //to here when finding only eigenvalues.
                end; 
                D[L]:=D[L] - P;
                E[L]:=G;
                E[M]:=0; 
                GoTo 1;
            end;
       end; 
    end;
end;