MODEが1の結果は、負の固有値によるものですから無視して下さい。計算は固有値の数を30出すように設定していますが、ERROR-IN-NORMからMODE=2からMODE=17までが信頼できる値となっています。上の表ではMODE=19まで表示しています。表を見ると、同じ固有振動数が2つづつ出ています。これは、x-yモードとx-zモードが同じだからです。MODE=10とMODE=13はねじれのモードだと思います。

計算精度を検証してみましょう。上の表を見ると、1次元の解と比較するとFEMの計算結果は、2桁まで合っています。１次元が正解とは言えませんが、FEMは良い結果を出していると思います。下にFEMの結果を1次元の解と比較してみました。縦軸が%で表した差((１次元-FEM)/１次元*100)で、横軸がX軸の要素数です。

X軸を100要素で分割した場合差は約7%、500要素で分割した場合差は約2%以内になっています。また、MODEが高次になるほど、MODE11を境に精度が悪化しています。これは、部分空間の[T]マトリクスの大きさを小さく(この場合30)とったことが原因だと考えられます。
下は固有ベクトルのプロットです。1次と2次モード(MODE02とMODE04)をプロットしてみました。予想した通りの結果になっています。