diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..75d4bbc 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,10 +1,40 @@
 #include "cotmatrix.h"
+#include <igl/doublearea.h>
+#include <cmath>
 
 void cotmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::SparseMatrix<double> & L)
 {
-  // Add your code here
+    L.resize(F.maxCoeff() + 1, F.maxCoeff() + 1);
+
+    Eigen::VectorXd area;
+    igl::doublearea(l, area);
+
+    std::vector<Eigen::Triplet<double>> tripletList;
+    tripletList.reserve(F.rows() * F.cols());
+
+    for (int i = 0; i < F.rows(); ++i) {
+        const Eigen::RowVector3d edge = l.row(i);
+
+        const double cot0 = (edge(1) * edge(1) + edge(2) * edge(2) - edge(0) * edge(0)) / (8 * area(i));
+        const double cot1 = (edge(0) * edge(0) + edge(2) * edge(2) - edge(1) * edge(1)) / (8 * area(i));
+        const double cot2 = (edge(1) * edge(1) + edge(0) * edge(0) - edge(2) * edge(2)) / (8 * area(i));
+
+        tripletList.push_back(Eigen::Triplet<double>(F(i, 0), F(i, 1), cot2));
+        tripletList.push_back(Eigen::Triplet<double>(F(i, 1), F(i, 2), cot0));
+        tripletList.push_back(Eigen::Triplet<double>(F(i, 2), F(i, 0), cot1));
+        tripletList.push_back(Eigen::Triplet<double>(F(i, 1), F(i, 0), cot2));
+        tripletList.push_back(Eigen::Triplet<double>(F(i, 2), F(i, 1), cot0));
+        tripletList.push_back(Eigen::Triplet<double>(F(i, 0), F(i, 2), cot1));
+    }
+
+    L.setFromTriplets(tripletList.begin(), tripletList.end());
+
+    for (int i = 0; i < L.rows(); ++i) {
+        // See issue #15
+        L.coeffRef(i, i) -= L.row(i).sum();
+    }
 }
 
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..06832ad 100644
--- a/src/massmatrix.cpp
+++ b/src/massmatrix.cpp
@@ -1,10 +1,26 @@
 #include "massmatrix.h"
+#include <igl/doublearea.h>
 
 void massmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
 {
-  // Add your code here
+    M.resize(F.maxCoeff() + 1);
+
+    Eigen::VectorXd area;
+    igl::doublearea(l, area);
+
+    for (int i = 0; i < M.rows(); i++) {
+      double sum = 0;
+
+      for (int j = 0; j < F.rows(); j++) {
+          if (F(j, 0) == i || F(j,1) == i || F(j, 2) == i) {
+              sum += area(j);
+          }
+      }
+
+      M.diagonal()[i] = sum / 3;
+    }
 }
 
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..a4d0966 100644
--- a/src/smooth.cpp
+++ b/src/smooth.cpp
@@ -1,4 +1,7 @@
 #include "smooth.h"
+#include "massmatrix.h"
+#include "cotmatrix.h"
+#include <igl/edge_lengths.h>
 
 void smooth(
     const Eigen::MatrixXd & V,
@@ -7,6 +10,21 @@ void smooth(
     double lambda,
     Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  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);
+
+    Eigen::SparseMatrix<double> A = -lambda * L;
+    for (int i = 0; i < M.rows(); ++i) {
+        A.coeffRef(i, i) += M.diagonal()(i);
+    }
+
+    Eigen::MatrixXd B = M * G;
+    Eigen::BiCGSTAB<Eigen::SparseMatrix<double>> solver(A);
+    U = solver.solve(B);
 }