The following code implements Simultaneous iteration for a matrix A∈(ℝ^(n×n)) on a list of vectors passed to the function as columns of a Tall Skinny matrix Q∈(ℝ^(n×m)).  The code checks for invariance of the subspace using theorem 6.1.3 and computes the eigenvalues using exercise 6.1.4 from our book i.e. it constructs the matrix B∈(ℝ^(m×m)) and check that
	AQ≈QB
using ||AQ-QB||<||A|| (10^(-7)). 
Once the subspace has converged it reports the eigenvalues of the small matrix B.

{n,m}={14,3};
A=RandomReal[{-1,1},{n,n}];
Q=RandomReal[{-1,1},{n,m}];
Q=Transpose[Orthogonalize[Transpose[Q]]];

ListLogPlot[
Table[
AQ=A.Q;
B=Transpose[Q].AQ;
err=Norm[AQ-Q.B];
Q=Transpose[Orthogonalize[Transpose[AQ]]];
err,
{250}],
PlotLabel->Eigenvalues[A,3]
]
MatrixForm[B]

Clear[SimIter]
SimIter[A_,QIn_]:= Module[
{err=∞,Q=QIn,n=0, AQ, B,MaxIter=300,Tol=Norm[A]*(10^(-7))},
err=∞;
Q=QIn;
While[
And[n++<MaxIter, norm>Tol],
AQ=A.Q;
Q=Transpose[Orthogonalize[Transpose[AQ]]];
B=Table[AQ[[i]].Q[[j]],{i,Qm},{j,Qm}];
norm=Norm[B.Q-AQ]; Print[norm];
Print[{MatrixForm[AQ], MatrixForm[B.Q]}];
];
{n,Q, B}
]