６元一次方程式の行列の条件数を求める方法

変数がθ２，θ３の２つがあるとして６元一次方程式を解く際にすべてのθ２，θ３の総当たりで解きたいのですが，方法はありますでしょうか．

θ２，θ３はそれぞれ‐９０～９０，‐９０～４５の範囲です．

eq1=a*sind(θ２+θ３)+b*cosd(θ２+θ３)+ c*sind(θ２)sind(θ２+θ３)+ d*sind(θ２)cosd(θ２+θ３)+ e*sind(θ２)cosd(θ２)+ f==I(91+θ２,91+θ３)

[A,B] = equationsToMatrix([eqn1,~eqn6], [a,~f]);

上記のような式が６つあり，Iは181行136列のデータ群です．

式のθ2，θ3を総当たりで計算して，各A行列の条件数Cond()最も低くなる6つのθ2，θ3の組み合わせを調べたいです．

現在ではfor文を用いて下記のように行おうとしていますがうまくいきません．

できれば，連立方程式を解く際に使用した角度6つがわかるようにしたいです．

for a3=-90:1:45
    for a2=90:-1:-90
                
         for b3=-90:1:45
            for b2=90:-1:-90
                
                for c3=-90:1:45
                    for c2=90:-1:90
                    
                        for d3=-90:1:45
                            for d2=90:-1:-90
                                
                                for e3=-90:1:45
                                    for e2=90:-1:-90
                                        
                                        for f3=-90:1:45
                                            for f2=90:-1:90
    
                                                syms a b c d e f

                                                     eqn1 =              a*sind(a3+a2).^2 ... 
                                                                        +b*sind(a2).^2 ...
                                                                        +c*sind(a2)*sind(a2+a3) ...
                                                                        +d*sind(a2+a3) ...
                                                                        +e*sind(a2) ...
                                                                        +f == I1_sim(91-a2,91+a3);

                                                    eqn2 =              a*sind(b2+b3).^2 ... 
                                                                        +b*sind(b2).^2 ...
                                                                        +c*sind(b2)*sind(b2+b3) ...
                                                                        +d*sind(b2+b3) ...
                                                                        +e*sind(b2) ...
                                                                        +f == I1_sim(91-b2,91+b3);

                                                    eqn3 =              a*sind(c2+c3).^2 ... 
                                                                        +b*sind(c2).^2 ...
                                                                        +c*sind(c2)*sind(c2+c3) ...
                                                                        +d*sind(c2+c3) ...
                                                                        +e*sind(c2) ...
                                                                        +f == I1_sim(91-c2,91+c3);

                                                    eqn4 =              a*sind(d2+d3).^2 ...
                                                                        +b*sind(d2).^2 ...
                                                                        +c*sind(d2)*sind(d2+d3) ...
                                                                        +d*sind(d2+d3) ...
                                                                        +e*sind(d2) ...
                                                                        +f == I1_sim(91-d2,91+d3);

                                                    eqn5 =              a*sind(e2+e3).^2 ... 
                                                                        +b*sind(e2).^2 ...
                                                                        +c*sind(e2)*sind(e2+e3) ...
                                                                        +d*sind(e2+e3) ...
                                                                        +e*sind(e2) ...
                                                                        +f == I1_sim(91-e2,91+e3);

                                                    eqn6 =              a*sind(f2+f3).^2 ... 
                                                                        +b*sind(f2).^2 ...
                                                                        +c*sind(f2)*sind(f2+f3) ...
                                                                        +d*sind(f2+f3) ...
                                                                        +e*sind(f2) ...
                                                                        +f == I1_sim(91-f2,91+f3);

                                                    [A,B] = equationsToMatrix([eqn1, eqn2, eqn3, eqn4, eqn5, eqn6], [a, b, c, d, e, f]);
                                                    X = linsolve(A,B);
                                                    Nolm(x,y,z)=cond(A);
                                                    x=x+1;
                                            end
                                            y=y+1;
                                        end
                                        z=z+1;
                                    end
                                    
                                end
                                z=z+1;

                            end
                            
                        end
                        z=z+1;
                    end
                    
                end
                z=z+1;
            end
            
        end
        z=z+1;
    end
    
end