diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..a8a51f5 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,4 +1,9 @@
 #include "cotmatrix.h"
+#include <vector>
+#include <cmath> // only for sqrt(.) and pow(.)
+
+// computes the cot of angle_ij given adj side lengths ki, kj, and the opposite side ij to angle_ij
+double compute_cot_from_lengths(double kj, double ki, double ij);
 
 void cotmatrix(
   const Eigen::MatrixXd & l,
@@ -6,5 +11,71 @@ void cotmatrix(
   Eigen::SparseMatrix<double> & L)
 {
   // Add your code here
+
+  std::vector<Eigen::Triplet<double> > L_vector;
+
+  int num_vertices = F.maxCoeff() + 1;
+
+  // duplicates get summed up when setting from triplets
+  // so, no need to interate on all possible edges
+
+  // iterate on the face matrix F to get all the existing edges
+  for (int f = 0; f < F.rows(); f++) {
+
+
+    // Vertex - index: i - 0; j - 1; k - 2
+
+    double len_ij = l(f, 2);
+    double len_ki = l(f, 1);
+    double len_kj = l(f, 0);
+
+    // handle the case where the edge e_ij exists
+    double cot_ij = compute_cot_from_lengths(len_kj, len_ki, len_ij);
+    double cot_kj = compute_cot_from_lengths(len_ki, len_ij, len_kj);
+    double cot_ki = compute_cot_from_lengths(len_ij, len_kj, len_ki);
+
+    L_vector.push_back(Eigen::Triplet<double>(F(f, 0), F(f, 1), 0.5*cot_ij));
+    L_vector.push_back(Eigen::Triplet<double>(F(f, 1), F(f, 0), 0.5*cot_ij));
+
+    L_vector.push_back(Eigen::Triplet<double>(F(f, 1), F(f, 2), 0.5*cot_kj));
+    L_vector.push_back(Eigen::Triplet<double>(F(f, 2), F(f, 1), 0.5*cot_kj));
+
+    L_vector.push_back(Eigen::Triplet<double>(F(f, 0), F(f, 2), 0.5*cot_ki));
+    L_vector.push_back(Eigen::Triplet<double>(F(f, 2), F(f, 0), 0.5*cot_ki));
+
+    // handle the case where i=j and edge exists 
+    // need negative of sum of all L_ik (for k != i) when edge exists between k and i
+    // exploting the property that duplicate entry gets added finally
+    L_vector.push_back(Eigen::Triplet<double>(F(f, 0), F(f, 0), -1.0*(0.5*cot_ij + 0.5*cot_ki)));
+    L_vector.push_back(Eigen::Triplet<double>(F(f, 1), F(f, 1), -1.0*(0.5*cot_ij + 0.5*cot_kj)));
+    L_vector.push_back(Eigen::Triplet<double>(F(f, 2), F(f, 2), -1.0*(0.5*cot_kj + 0.5*cot_ki)));
+
+    // rest of the cases, it's zero. so nothing to do.
+    
+  }
+
+  L.resize(num_vertices, num_vertices);
+  L.setFromTriplets(L_vector.begin(), L_vector.end());
+
+}
+
+// computes the cot of angle_ij given adj side lengths ki, kj, and the opposite side ij to angle_ij
+double compute_cot_from_lengths(double kj, double ki, double ij) {
+
+  // first compute the area using Heron's formula
+  double semi_perimeter = (kj + ki + ij)/2.0;
+  double area_sq = semi_perimeter * (semi_perimeter - kj) * (semi_perimeter - ki) * (semi_perimeter - ij) * 1.0;
+  double area = sqrt(area_sq);
+
+  // use area = 0.5 ab Sin(C) to find Sine angles
+  double sine_ij = (2.0*area)/(kj*ki*1.0);
+
+  // use cosine rule to compute cos(angle_ij)
+  double cosine_ij = (kj*kj + ki*ki - ij*ij)/(2.0*kj*ki);
+
+  // cot_ij
+  double cot_ij = cosine_ij/sine_ij;
+
+  return cot_ij;
 }
 
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..54c8815 100644
--- a/src/massmatrix.cpp
+++ b/src/massmatrix.cpp
@@ -1,4 +1,7 @@
 #include "massmatrix.h"
+#include <cmath>
+
+double compute_area_herons_formula(double a, double b, double c);
 
 void massmatrix(
   const Eigen::MatrixXd & l,
@@ -6,5 +9,38 @@ void massmatrix(
   Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
 {
   // Add your code here
+  int num_vertices = F.maxCoeff() + 1;
+  
+  Eigen::VectorXd M_vec(num_vertices);
+  M_vec.setZero();
+
+  // iterate over the faces
+  for (int f = 0; f < F.rows(); f++) {
+
+    // get the lengths
+    double a = l(f, 0);
+    double b = l(f, 1);
+    double c = l(f, 2);
+
+    // get the area of the triangle
+    double area = compute_area_herons_formula(a, b, c);
+
+    M_vec(F(f, 0)) += 0.333 * area;
+    M_vec(F(f, 1)) += 0.333 * area;
+    M_vec(F(f, 2)) += 0.333 * area;
+  }
+
+  M.resize(num_vertices);
+  M = M_vec.asDiagonal();
+}
+
+
+double compute_area_herons_formula(double a, double b, double c) {
+
+  double semi_per = (a+b+c)/2.0;
+  double area_sq = semi_per * (semi_per - a) * (semi_per - b) * (semi_per - c);
+  double area = sqrt(area_sq);
+
+  return area;
 }
 
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..058ed40 100644
--- a/src/smooth.cpp
+++ b/src/smooth.cpp
@@ -1,4 +1,9 @@
 #include "smooth.h"
+#include "cotmatrix.h"
+#include "massmatrix.h"
+#include <igl/edge_lengths.h>
+#include <Eigen/SparseCholesky>
+#include <vector>
 
 void smooth(
     const Eigen::MatrixXd & V,
@@ -8,5 +13,35 @@ void smooth(
     Eigen::MatrixXd & U)
 {
   // Replace with your code
-  U = G;
+  // U = G;
+
+  Eigen::MatrixXd l;
+  igl::edge_lengths(V, F, l);
+
+  Eigen::SparseMatrix<double> L;
+  cotmatrix(l, F, L);
+
+  Eigen::DiagonalMatrix<double,Eigen::Dynamic> M;
+  massmatrix(l, F, M);
+
+  // create a sparse version of M as DiagonalMatrix is not able to add or type cast to a SparseMatrix
+  Eigen::VectorXd M_vec = M.diagonal();
+  Eigen::SparseMatrix<double> M_sparse;
+  std::vector<Eigen::Triplet<double> > M_sparse_vector;
+  for (int k = 0; k<M_vec.size(); k++) {
+    M_sparse_vector.push_back(Eigen::Triplet<double>(k, k, M_vec(k)));
+  }
+  M_sparse.resize(M_vec.size(), M_vec.size());
+  M_sparse.setFromTriplets(M_sparse_vector.begin(), M_sparse_vector.end());
+  
+
+  U.resize(G.rows(), G.cols());
+
+  // Cholesky decomposition and solving
+  // The usual SimplicalLLT is very very slow. Resorting to LDLT version of Cholesky
+  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
+
+  solver.compute(M_sparse - lambda * L);
+  U = solver.solve(M_sparse * G);
+
 }