diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp
index a87325e..9369245 100644
--- a/include/tug/Grid.hpp
+++ b/include/tug/Grid.hpp
@@ -132,6 +132,13 @@ class Grid {
          */
         void setDomain(int domainRow, int domainCol);
 
+        /**
+         * @brief Gets the delta value for 1D-Grid. Grid must be one dimensional.
+         * 
+         * @return double Delta value. 
+         */
+        double getDelta();
+
         /**
          * @brief Gets the delta value in x-direction. 
          * 
diff --git a/src/BTCSv2.cpp b/src/BTCSv2.cpp
index 3fc5c69..8b7ed21 100644
--- a/src/BTCSv2.cpp
+++ b/src/BTCSv2.cpp
@@ -1,9 +1,17 @@
+/**
+ * @file BTCS.cpp
+ * @brief Implementation of heterogenous BTCS (backward time-centered space) solution
+ *        of diffusion equation in 1D and 2D space.
+ * 
+ */
+
 #include "FTCS.cpp"
 #include <tug/Boundary.hpp>
 
 using namespace Eigen;
 
 
+// calculates coefficient for left boundary in constant case
 static tuple<double, double> calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int &rowIndex, double &sx) {
     double centerCoeff;
     double rightCoeff;
@@ -16,6 +24,7 @@ static tuple<double, double> calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int
 }
 
 
+// calculates coefficient for left boundary in closed case
 static tuple<double, double> calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int &rowIndex, double &sx) {
     double centerCoeff;
     double rightCoeff;
@@ -27,6 +36,7 @@ static tuple<double, double> calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int &r
 }
 
 
+// calculates coefficient for right boundary in constant case
 static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha, int &rowIndex, int &n, double &sx) {
     double leftCoeff;
     double centerCoeff;
@@ -39,6 +49,7 @@ static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha, int
 }
 
 
+// calculates coefficient for right boundary in closed case
 static tuple<double, double> calcRightBoundaryCoeffClosed(MatrixXd &alpha, int &rowIndex, int &n, double &sx) {
     double leftCoeff;
     double centerCoeff;
@@ -50,6 +61,7 @@ static tuple<double, double> calcRightBoundaryCoeffClosed(MatrixXd &alpha, int &
 }
 
 
+// creates coefficient matrix for next time step from alphas in x-direction
 static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, vector<BoundaryElement> bcLeft, vector<BoundaryElement> bcRight, int numCols, int rowIndex, double sx) {
 
     // square matrix of column^2 dimension for the coefficients
@@ -102,6 +114,7 @@ static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, vector<BoundaryEl
 }
 
 
+// calculates explicity concentration at top boundary in constant case
 static double calcExplicitConcentrationsTopBoundaryConstant(MatrixXd &concentrations, 
                             MatrixXd &alpha, vector<BoundaryElement> &bcTop, int &rowIndex, int &i, double &sy) {
     double c;
@@ -120,6 +133,7 @@ static double calcExplicitConcentrationsTopBoundaryConstant(MatrixXd &concentrat
 }
 
 
+// calculates explicit concentration at top boundary in closed case
 static double calcExplicitConcentrationsTopBoundaryClosed(MatrixXd &concentrations, 
                             MatrixXd &alpha, int &rowIndex, int &i, double &sy) {
     double c;
@@ -136,6 +150,7 @@ static double calcExplicitConcentrationsTopBoundaryClosed(MatrixXd &concentratio
 }
 
 
+// calculates explicit concentration at bottom boundary in constant case
 static double calcExplicitConcentrationsBottomBoundaryConstant(MatrixXd &concentrations, 
                             MatrixXd &alpha, vector<BoundaryElement> &bcBottom, int &rowIndex, int &i, double &sy) {
     double c;
@@ -154,6 +169,7 @@ static double calcExplicitConcentrationsBottomBoundaryConstant(MatrixXd &concent
 }
 
 
+// calculates explicit concentration at bottom boundary in closed case
 static double calcExplicitConcentrationsBottomBoundaryClosed(MatrixXd &concentrations, 
                             MatrixXd &alpha, int &rowIndex, int &i, double &sy) {
     double c;
@@ -170,6 +186,7 @@ static double calcExplicitConcentrationsBottomBoundaryClosed(MatrixXd &concentra
 }
 
 
+// creates a solution vector for next time step from the current state of concentrations
 static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY, 
                                         vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight, 
                                         vector<BoundaryElement> &bcTop, vector<BoundaryElement> &bcBottom, 
@@ -238,6 +255,8 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX,
 }
 
 
+// solver for linear equation system; A corresponds to coefficient matrix,
+// b to the solution vector
 static VectorXd solve(SparseMatrix<double> &A, VectorXd &b) {
     SparseLU<SparseMatrix<double>> solver;
     solver.analyzePattern(A);
@@ -250,7 +269,7 @@ static VectorXd solve(SparseMatrix<double> &A, VectorXd &b) {
 // BTCS solution for 1D grid 
 static void BTCS_1D(Grid &grid, Boundary &bc, double &timestep) {
     int length = grid.getLength();
-    double sx = timestep / (grid.getDeltaCol() * grid.getDeltaCol()); // TODO create method getDelta() for 1D case
+    double sx = timestep / (grid.getDelta() * grid.getDelta());
 
     VectorXd concentrations_t1(length);
 
@@ -312,7 +331,7 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
         row_t1 = solve(A, b);
         
         for (int j = 0; j < colMax; j++) {
-            concentrations_t1(i,j) = row_t1(j); // TODO Eigen seq ??
+            concentrations_t1(i,j) = row_t1(j); // can potentially be improved by using Eigen method
         }
         
     }
@@ -329,7 +348,7 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
         row_t1 = solve(A, b);
 
         for (int j = 0; j < rowMax; j++) {
-            concentrations(i,j) = row_t1(j); // TODO Eigen seq ??
+            concentrations(i,j) = row_t1(j); // can potentially be improved by using Eigen method
         }
     }
     concentrations.transposeInPlace();
diff --git a/src/Grid.cpp b/src/Grid.cpp
index 3327fb9..66e95ee 100644
--- a/src/Grid.cpp
+++ b/src/Grid.cpp
@@ -145,6 +145,14 @@ void Grid::setDomain(int domainRow, int domainCol) {
     this->deltaCol = double(this->domainCol)/double(this->col);
 }
 
+double Grid::getDelta() {
+    if (dim != 1) {
+        throw_invalid_argument("Grid is not one dimensional, you should probably use the 2D delta getters");
+    }
+
+    return this->deltaCol;
+}
+
 double Grid::getDeltaCol() {
     return this->deltaCol;
 }
diff --git a/src/README.md b/src/README.md
new file mode 100644
index 0000000..420608f
--- /dev/null
+++ b/src/README.md
@@ -0,0 +1,18 @@
+<h1>src Directory</h1>
+
+<pre>
+src/  
+    ├── CMakeFiles/
+    ├── Boundary.cpp  
+    ├── BoundaryCondition.cpp 
+    ├── BTCS.cpp  
+    ├── BTCSv2.cpp
+    ├── CMakeLists.txt  
+    ├── FTCS.cp
+    ├── Grid.cpp 
+    ├── README.md
+    ├── Simulation.cpp
+    ├── Solver.cpp  
+    └── TugUtils.hpp 
+
+</pre>
\ No newline at end of file