Решение СЛАУ методом Гаусса

//---------------------------------------------------------------------------
#pragma hdrstop
#include "Unit2.h"
#include <math.h>
//---------------------------------------------------------------------------
#pragma package(smart_init)
class gauss
{
    double** matr;
    double* vek;
    double* otv;
    double* mas;
    int n;
    void func();
public:
    gauss(double** c,double* d,int N);
    friend double* func1(gauss m);
};
gauss::gauss(double** c,double* d,int N)
{
    n = N;
    matr = new double*[n];
    for (int i = 0; i < n; i++)
    {
 
        matr[i] = new double[n];
        for (int j = 0; j < n; j++)
        {
            matr[i][j] = c[i][j];
        }
    }
    vek = new double[n];
    for (int i = 0; i < n; i++)
    {
        vek[i] = d[i];
    }
    mas = new double[n];
    otv = new double[n];
    func();
}
void gauss::func()
{
int l, p, M;
double a, b, c;
// присвоение компонентам массива перестановок mas(k) исходных значений
       for (int k=0; k<n; k++)
       {
          mas[k] = k;
       }
//Hахождение индекса p
for (int k=0; k<n; k++)
 {
 
        a=0;
          for (int i=k; i<n; i++)
          {
             l=mas[i];
           if ( fabs (matr[l][k]) < a )
             continue;
             M = l;
             p = i;
             a = fabs ( matr[l][k]  );
          }
      b = matr[M][k]; //Замена местами значения mas[p]  и mas[k], если mas[p] не равно mas[k]
   if (k!=p)
   {
    mas[p] = mas[k];
    mas[k] = M;
   };
   if (b == 0)
   {
    break;
   };
        for (int j = k; j<n; j++)    //Исключение переменной otv[k] (прямой ход метода Гаусса)
        {
          matr[M][j] = matr[M][j] / b;
        }
        vek[M] = vek[M]/b;
 
        for (int i = k+1; i<n; i++)
        {
         l = mas[i];
         c = matr[l][k];
              for (int j = k; j<n; j++)
              {
                matr[l][j] = matr[l][j] - c * matr[M][j];
              }
          vek[l] = vek[l] - c * vek[M]; 
         }
   if ( matr[l][n-1] == 0 ) 
   {
    break;
   };
 }
/*Обратный ход метода Гаусса*/
double sum;
       for ( int k = n-1; k>=0; k--)
       {
          l = mas[k];
          sum = 0;
             for (int j = k+1; j<n; j++)
             {
                sum = sum + matr[l][j] * otv[j];
             }
         otv[k] = vek[l] - sum;
        }
 
}
double* func1(gauss m)
{
    return m.otv;
}