Program to solve a system of linear equations in C++
<p>I am testing this code for solving linear systems with this simple 2-equation system (in matrix form "Mat[2][3]"), but when I execute it, I obtain the following result, which does not agree with the coefficients I have introduced in the system Matrix:</p>

<p>CODE:</p>

<pre><code>//Gauss Elimination
#include &lt;iostream&gt;
#include &lt;cmath&gt;
#include &lt;vector&gt;
#include &lt;fstream&gt;
#include &lt;iterator&gt;
#include &lt;string&gt;
#include &lt;iomanip&gt;
using namespace std;

int main() {

    double Mat[2][3];
    vector&lt;double&gt; q(2);
    int Nx = 2, ii, jj, kk;

    Mat[1][1] = 2.0;
    Mat[1][2] = 3.0;
    Mat[1][3] = 2.5;
    Mat[2][1] = 1.3;
    Mat[2][2] = 3.0;
    Mat[2][3] = 2.5;         

    cout &lt;&lt; "Matrix: " &lt;&lt; endl;     
    for (ii = 0; ii &lt; Nx; ii++) {
        for (jj = 0; jj &lt; Nx + 1; jj++) {
            cout &lt;&lt; Mat[ii][jj] &lt;&lt; "   ";
            cout &lt;&lt; endl;
        }
    }
    // Triangularization
    for (ii = 0; ii &lt; Nx - 1; ii++)
        for (kk = ii + 1; kk &lt; Nx; kk++)
        {
            double t = Mat[kk][ii] / Mat[ii][ii];
            for (jj = 0; jj &lt;= Nx; jj++)
                Mat[kk][jj] = Mat[kk][jj] - t * Mat[ii][jj];
        }   // Resolution
    for (ii = Nx - 1; ii &gt;= 0; ii--)
    {
        q[ii] = Mat[ii][Nx];
        for (jj = Nx - 1; jj &gt; ii; jj--)                        
            q[ii] = q[ii] - Mat[ii][jj] * q[jj];
        q[ii] = q[ii] / Mat[ii][ii];
    }

    cout &lt;&lt; "Solution of the system: " &lt;&lt; endl;
    cout &lt;&lt; q[1] &lt;&lt; endl;
    cout &lt;&lt; q[2] &lt;&lt; endl; ////////////////////////////////////////////////////////////////////////


    return 0;
}
</code></pre>

<p>RESULT:</p>

<pre><code>Matrix: 
0   2.07496e-317   6.95314e-310   
0   2   3   
Solution of the system: 
-nan
0
</code></pre>