Matrix SqMatrix::solve(Matrix b) {
  assert(b.columns()==1);
  int dim=this->dimension();

  Matrix y (dim, 1);
  Matrix x (dim, 1);

  //forward substitution algorithm - solve for y, Ly=b

  for(int i=0; i<dim; i++) {
    double s=0;  //the sum in the equation
    for(int j=0; j<i; j++) {
      s+=(*this)(i, j) * y(j,0);
    }
    y(i, 0)=(b(i,0)-s);  // / (*this)(i, i); das sollte gleich eins werden bei L
    if(y(i,0)==0.0) y(i,0)=0;
  }

  //backward substitution - Ux=y, solve for x

  for(int i=this->dimension()-1; i>=0; i--) {
    double s=0;
    for(int j=i; j<this->dimension()-1; j++) {
      s+=(*this)(i, j)*x(j,0);
    }
    x(i,0)=(b(i,0)-s) / (*this)(i, i);
    if(x(i,0)==0.0) x(i,0)=0;
  }

  return x;


}