diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp
index 1319677..6a6f973 100644
--- a/include/tug/Grid.hpp
+++ b/include/tug/Grid.hpp
@@ -43,7 +43,7 @@ public:
    * must have correct dimensions as defined in row and col. (Or length, in 1D
    * case).
    */
-  void setConcentrations(Eigen::MatrixXd concentrations);
+  void setConcentrations(const Eigen::MatrixXd &concentrations);
 
   /**
    * @brief Gets the concentrations matrix for a Grid.
@@ -51,7 +51,7 @@ public:
    * @return MatrixXd An Eigen3 matrix holding the concentrations and having the
    *                  same dimensions as the grid.
    */
-  const Eigen::MatrixXd getConcentrations();
+  const Eigen::MatrixXd &getConcentrations() { return this->concentrations; }
 
   /**
    * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one
@@ -60,7 +60,7 @@ public:
    * @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients.
    * Matrix columns must have same size as length of grid.
    */
-  void setAlpha(Eigen::MatrixXd alpha);
+  void setAlpha(const Eigen::MatrixXd &alpha);
 
   /**
    * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two
@@ -71,7 +71,7 @@ public:
    * @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in
    * y-direction. Matrix must be of same size as the grid.
    */
-  void setAlpha(Eigen::MatrixXd alphaX, Eigen::MatrixXd alphaY);
+  void setAlpha(const Eigen::MatrixXd &alphaX, const Eigen::MatrixXd &alphaY);
 
   /**
    * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one
@@ -79,7 +79,7 @@ public:
    *
    * @return MatrixXd A matrix with 1 row holding the alpha coefficients.
    */
-  const Eigen::MatrixXd getAlpha();
+  const Eigen::MatrixXd &getAlpha();
 
   /**
    * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid.
@@ -87,7 +87,7 @@ public:
    *
    * @return MatrixXd A matrix holding the alpha coefficients in x-direction.
    */
-  const Eigen::MatrixXd getAlphaX();
+  const Eigen::MatrixXd &getAlphaX();
 
   /**
    * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid.
@@ -95,7 +95,7 @@ public:
    *
    * @return MatrixXd A matrix holding the alpha coefficients in y-direction.
    */
-  const Eigen::MatrixXd getAlphaY();
+  const Eigen::MatrixXd &getAlphaY();
 
   /**
    * @brief Gets the dimensions of the grid.
@@ -165,12 +165,12 @@ public:
 
 private:
   int col;                        // number of grid columns
-  int row;                        // number of grid rows
+  int row{1};                     // number of grid rows
   int dim;                        // 1D or 2D
   double domainCol;               // number of domain columns
-  double domainRow;               // number of domain rows
+  double domainRow{0};            // number of domain rows
   double deltaCol;                // delta in x-direction (between columns)
-  double deltaRow;                // delta in y-direction (between rows)
+  double deltaRow{0};             // delta in y-direction (between rows)
   Eigen::MatrixXd concentrations; // Matrix holding grid concentrations
   Eigen::MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction
   Eigen::MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction
diff --git a/src/Grid.cpp b/src/Grid.cpp
index dad79e9..c6573e4 100644
--- a/src/Grid.cpp
+++ b/src/Grid.cpp
@@ -3,42 +3,38 @@
 #include <iostream>
 #include <tug/Grid.hpp>
 
-Grid::Grid(int length) {
+constexpr double MAT_INIT_VAL = 0;
+
+Grid::Grid(int length) : col(length), domainCol(length) {
   if (length <= 3) {
     throw_invalid_argument(
         "Given grid length too small. Must be greater than 3.");
   }
 
-  this->row = 1;
-  this->col = length;
-  this->domainCol = length; // default: same size as length
-  this->deltaCol = double(this->domainCol) / double(this->col); // -> 1
   this->dim = 1;
+  this->deltaCol = double(this->domainCol) / double(this->col); // -> 1
 
-  this->concentrations = Eigen::MatrixXd::Constant(1, col, 20);
-  this->alphaX = Eigen::MatrixXd::Constant(1, col, 1);
+  this->concentrations = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL);
+  this->alphaX = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL);
 }
 
-Grid::Grid(int row, int col) {
+Grid::Grid(int _row, int _col)
+    : row(_row), col(_col), domainRow(_row), domainCol(_col) {
   if (row <= 3 || col <= 3) {
     throw_invalid_argument(
         "Given grid dimensions too small. Must each be greater than 3.");
   }
 
-  this->row = row;
-  this->col = col;
-  this->domainRow = row; // default: same size as row
-  this->domainCol = col; // default: same size as col
+  this->dim = 2;
   this->deltaRow = double(this->domainRow) / double(this->row); // -> 1
   this->deltaCol = double(this->domainCol) / double(this->col); // -> 1
-  this->dim = 2;
 
-  this->concentrations = Eigen::MatrixXd::Constant(row, col, 20);
-  this->alphaX = Eigen::MatrixXd::Constant(row, col, 1);
-  this->alphaY = Eigen::MatrixXd::Constant(row, col, 1);
+  this->concentrations = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL);
+  this->alphaX = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL);
+  this->alphaY = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL);
 }
 
-void Grid::setConcentrations(Eigen::MatrixXd concentrations) {
+void Grid::setConcentrations(const Eigen::MatrixXd &concentrations) {
   if (concentrations.rows() != this->row ||
       concentrations.cols() != this->col) {
     throw_invalid_argument(
@@ -48,9 +44,7 @@ void Grid::setConcentrations(Eigen::MatrixXd concentrations) {
   this->concentrations = concentrations;
 }
 
-const Eigen::MatrixXd Grid::getConcentrations() { return this->concentrations; }
-
-void Grid::setAlpha(Eigen::MatrixXd alpha) {
+void Grid::setAlpha(const Eigen::MatrixXd &alpha) {
   if (dim != 1) {
     throw_invalid_argument("Grid is not one dimensional, you should probably "
                            "use 2D setter function!");
@@ -63,7 +57,8 @@ void Grid::setAlpha(Eigen::MatrixXd alpha) {
   this->alphaX = alpha;
 }
 
-void Grid::setAlpha(Eigen::MatrixXd alphaX, Eigen::MatrixXd alphaY) {
+void Grid::setAlpha(const Eigen::MatrixXd &alphaX,
+                    const Eigen::MatrixXd &alphaY) {
   if (dim != 2) {
     throw_invalid_argument("Grid is not two dimensional, you should probably "
                            "use 1D setter function!");
@@ -81,7 +76,7 @@ void Grid::setAlpha(Eigen::MatrixXd alphaX, Eigen::MatrixXd alphaY) {
   this->alphaY = alphaY;
 }
 
-const Eigen::MatrixXd Grid::getAlpha() {
+const Eigen::MatrixXd &Grid::getAlpha() {
   if (dim != 1) {
     throw_invalid_argument("Grid is not one dimensional, you should probably "
                            "use either getAlphaX() or getAlphaY()!");
@@ -90,7 +85,7 @@ const Eigen::MatrixXd Grid::getAlpha() {
   return this->alphaX;
 }
 
-const Eigen::MatrixXd Grid::getAlphaX() {
+const Eigen::MatrixXd &Grid::getAlphaX() {
   if (dim != 2) {
     throw_invalid_argument(
         "Grid is not two dimensional, you should probably use getAlpha()!");
@@ -99,7 +94,7 @@ const Eigen::MatrixXd Grid::getAlphaX() {
   return this->alphaX;
 }
 
-const Eigen::MatrixXd Grid::getAlphaY() {
+const Eigen::MatrixXd &Grid::getAlphaY() {
   if (dim != 2) {
     throw_invalid_argument(
         "Grid is not two dimensional, you should probably use getAlpha()!");