Метод Гаусса на с++

Нашёл этот алгоритм для решения линейных уравнений через метод Гаусса

#include "stdafx.h"
#include <iostream>
#include <fstream>
using namespace std;
int n, i, j, k;
double d, s;
int main()
{
    ifstream fin("input.txt");
    cout << "vvedite N: " << endl;
    fin >> n;
    double **a = new double *[n];
    for (i = 0; i <= n; i++)
        a[i] = new double[n];
    double **a1 = new double *[n];
    for (i = 0; i <= n; i++)
        a1[i] = new double[n];
    double *b = new double[n];
    double *x = new double[n];
    for (i = 1; i <= n; i++)
    {
        for (j = 1; j <= n; j++)
        {
            fin >> a[i][j];
            a1[i][j] = a[i][j];
            cout << "a[ " << i << "," << j << "]= "<<a[i][j]<<endl;
        }
    }
    for (i = 1; i <= n; i++)
    {
        fin >> b[i];
        cout << "b,[ " << i << "]= "<<b[i]<<endl;
    }
    for (k = 1; k <= n; k++) 
    {
        for (j = k + 1; j <= n; j++)
        {
            d = a[j][k] / a[k][k]; 
            for (i = k; i <= n; i++)
            {
                a[j][i] = a[j][i] - d * a[k][i]; 
            }
            b[j] = b[j] - d * b[k]; 
        }
    }
    for (k = n-2; k >= 1; k--) 
    {
        d = 0;
        for (j = k + 1; j <= n; j++)
        {
            s = a[k][j] * x[j]; 
            d = d + s; 
        }
        x[k] = (b[k] - d) / a[k][k]; 
    }
    cout << "Korni sistemy: " << endl;
    ofstream out;
    out.open("output.txt");
    for (i = 1; i <= n; i++)
    {
        out <<'x'<<i<<'='<< x[i] << endl;
        cout << 'x' << i << '=' << x[i] << endl;
    }
    system("pause");
    return 0;
}

Но при вводе например такой матрицы

1 0 0 0 0 | 1

0 1 0 0 0 | 1

0 0 0 1 0 | 1

0 0 1 0 0 | 1

0 0 0 0 1 | 1

выдаётся -nan(ind). Видимо где-то происходит деление на ноль. Как решить проблему?