C++ (g++ 4.8.1) memory corruption && corrupted unsorted chunks

//Matrix.h
#define uint unsigned int                        // Писать каждый раз unsigned int нет желания
class matrix                                     // Класс для работы с матрицами
{
 private:                                        
         double**_array;                         // Массив массивов через указатели этакая матрешка, представляет нашу матрицу
         uint _width;                            // Размерности  
         uint _height;                           // нашей матрицы
  public:
         matrix(matrix&,int=0,int=0);            // Конструктор копии с возможностью сдвига (обрезки матрицы)
         matrix(uint,uint);                      // Инициализирующий конструктор
         ~matrix();                              // Деструктор где чистим память занимаемую _array

         double det();                           // Считает определитель
         double maxElementInRow(uint,uint);      // Находит максимальный элемент в строке
         double maxElementInColumn(uint,uint);   // Находит максимальный элемент в столбце
         double* operator[](uint);               // Индексация по матрице
         void changeRows(uint,uint);             // Поменять строки местами
         void changeColumns(uint,uint);          // Поменять столбцы местами
         uint getIndexOfRowEl(double,uint);      // Получение индекса элемента из указанного столбца
         uint getIndexOfColumnEl(double,uint);   // Получение индекса строки элемента по столбцу
         uint width();                           // Получение приватных полей _width _height
         uint height();                          
};

//Matrix.cpp
#include "Matrix.h"

matrix::matrix(matrix& matrix,int dWidth,int dHeight):
                                                       _array    (new double*[matrix._width+dWidth]),
                                                       _width    (matrix._width+dWidth),
                                                       _height   (matrix._height+dHeight)
{
    for(uint i=0;i<_width;i++)
    {  
        _array[i]=new double[_height];

        for(uint j=0;j<matrix._height;j++)
            {
                if(i<matrix._width)
                _array[i][j]=matrix[i][j];
            }
    }
}

matrix::matrix(uint width, uint height):
                                        _array    (new double*[width]),
                                        _width    (width),
                                        _height   (height)
{
    for(uint i=0;i<_width;i++)
        {
            _array[i]=new double[_height];
        }
}

 matrix::~matrix() // Освобождаем память
{
    for(uint i=0;i<_width;i++)
    {
        delete[] _array[i];
    }

    delete[] _array;
}

double* matrix::operator[](uint i)
{
    return _array[i];
}

uint matrix::width()
{
    return _width;
}

uint matrix::height()
{
    return _height;
}

void matrix::changeRows(uint firstRow,uint secondRow)
 {
    double temp;

    for(uint i=0;i<_width;i++)
    {
           temp=_array[i][secondRow];
           _array[i][secondRow]=_array[i][firstRow];
           _array[i][firstRow]=temp;
    }
 }

 void matrix::changeColumns(uint firstRow,uint secondRow)
 {
   double temp;

   for(uint i=0;i<_width;i++)
   {
           temp=_array[secondRow][i];
           _array[secondRow][i]=_array[firstRow][i];
           _array[firstRow][i]=temp;
   }
 }

 double matrix::maxElementInRow(uint row,uint offset)
 {
    double max=_array[offset][row];

    for(uint i=offset;i<_width-1;i++)
    {
        if(max<_array[i+1][row])
            max=_array[i+1][row];
    }

    return max;
 }

 uint matrix::getIndexOfRowEl(double element,uint row)
 {
    uint index=0;

    for(uint i=0;i<_width;i++)
    {
        if(_array[i][row]==element)
            index=i;
    }

    return index;
 }

 uint matrix::getIndexOfColumnEl(double element,uint column)
 {
    uint index=0;

    for(uint i=0;i<_height;i++)
    {
        if(_array[column][i]==element)
            index=i;
    }

    return index;
 }

 double matrix::maxElementInColumn(uint column,uint offset)
 {
    double max=_array[column][offset];

    for(uint i=offset;i<_width-1;i++)
    {
        if(max<_array[column][i+1])
            max=_array[column][i+1];
    }

    return max;
 }

double matrix::det()             // Перегоняем копию матрицы объекта в треугольный вид (Гауссом)  
{                                // и умножаем диагональные эл-ты что даёт на определитель
    matrix temp(*this);
    double det=1;

    for(uint i=0; i<temp._width-1; i++)
    {
        if(temp[i][i]!=temp.maxElementInColumn(i,i))
        {
            temp.changeRows(i,temp.getIndexOfColumnEl(temp.maxElementInColumn(i,i),i));
        }

        double firstVal=temp[i][i],
               SecondVal;

        for(uint j=i+1;j<temp._height;j++)
        {
            SecondVal=temp[i][j];

            for(uint k=0;k<temp._width;k++)
            {
                temp[k][j]=temp[k][j]-temp[k][i]*(SecondVal/firstVal);
            }
        }
    }

    for(uint i=0;i<temp._width;i++)
    {
        det*=temp[i][i];
    }

    return det;
}

//MathMethods.h
#include "Matrix.h"

class MathMethods
{  
    public:
            static double* GetSolutionViaGauss(matrix); // Решение СЛАУ методом Гаусса
            static double* Interpolation (double (*)(double),uint,double,double); // Получение интерполирующего полинома,
            //принимает указатель на функцию(что-то типа делегата если это можно им назвать)

};

//MathMethods.cpp
#include <stdlib.h>
#include <math.h>
#include "MathMethods.h"

double* MathMethods::GetSolutionViaGauss(matrix slau)
{
//  if(matrix(slau,-1).det()!=0)
        for(int i=0; i<slau.width()-2; i++)
        {
            if(slau[i][i]!=slau.maxElementInColumn(i,i))
                slau.changeRows(i,slau.getIndexOfColumnEl(slau.maxElementInColumn(i,i),i));

            double firstVal=slau[i][i],
                   SecondVal;

            for(int j=i+1;j<slau.height();j++)
           {
                SecondVal=slau[i][j];

                for(int k=0;k<slau.width();k++)
                    slau[k][j]=slau[k][j]-slau[k][i]*(SecondVal/firstVal);
           }
       }

        double* solution=new double (slau.height());

        for(int i=0;i<slau.height();i++)
            {
                solution[i]=0;
            }

        for(int i=slau.height()-1;i>-1;i--)
        {
         double temp=0;

         for(int j=0;j<slau.height();j++)
             {
                temp-=slau[j][i]*solution[j];
             }

         solution[i]=(temp+slau[slau.width()-1][i])/slau[i][i];
        }

        return solution;
    }

    double* MathMethods::Interpolation (double (*func)(double),uint n,double beginX,double endX)
    {
        uint N = n+2;
        matrix slau(N+1,N);
        double* arg = new double [N];

        arg[0] = beginX;
        arg[1] = endX;

        for(uint i=2; i<N; i++)
        {
            arg[i]=(double)rand()/(double)RAND_MAX*(endX-beginX)+beginX;
        }

        uint k = N;

        for(uint i=0; i<slau.width(); i++,k--)
        {
            for(uint j=0; j<slau.height(); j++)
            {
                if(i==N)
                {
                    slau[i][j]=1;
                }
                else if (i==N+1)
                {
                    slau[i][j] = func(arg[j]);
                }
                     else
                     {
                        slau[i][j]=pow(arg[j],k);
                     }
            }

        }
        delete[] arg;
        return MathMethods::GetSolutionViaGauss(slau);
    }

//Program.cpp с точкой входа
#include <iostream>
#include <string>
#include "MathMethods.h"

int main ()
{
    uint n;
    double beginX,endX;
    std::string result="F(x)=";
    std::cout<<"Введите количество узлов интерполяции:n N=";
    std::cin>>n;
    std::cout<<"nВведите отрезок функции для интерполирования:n От ";
    std::cin>>beginX;
    std::cout<<" До ";
    std::cin>>endX;
    std::cout<<"nИнтерполяционная функция:n";
    double* solution = MathMethods::Interpolation([](double x){return 2*x+5;},n,beginX,endX);
    n+=2;
    uint k=n-1;
   
    for(uint i=0; i<n; i++,k--)
    {
        if(solution[i]>0 && i>0)
        {
            result+="+"+std::to_string(solution[i])+"*X^"+std::to_string(k);
        }
        else
        {
            result+=std::to_string(solution[i])+"*X^"+std::to_string(k);
        }
    }
    std::cout<<result<<"nn>> Для выхода нажмите любую клавишу <<n";
    std::cin.get();
    return 0;
}