diff --git a/include/cotmatrix.h b/include/cotmatrix.h
index dd92cca..a6a2431 100644
--- a/include/cotmatrix.h
+++ b/include/cotmatrix.h
@@ -2,6 +2,7 @@
 #define COTMATRIX_H
 #include <Eigen/Core>
 #include <Eigen/Sparse>
+#include <cmath>
 // Construct the "cotangent Laplacian" for a mesh with edge lengths `l`. Each
 // entry in the output sparse, symmetric matrix `L` is given by:
 //
diff --git a/include/massmatrix.h b/include/massmatrix.h
index b822a5e..6d911d7 100644
--- a/include/massmatrix.h
+++ b/include/massmatrix.h
@@ -2,8 +2,9 @@
 #define MASSMATRIX_H
 #include <Eigen/Core>
 #include <Eigen/Sparse>
+#include <igl/doublearea.h>
 // Construct the diagonal(ized) mass matrix for a mesh with edge lengths `l`.
-// Each enetry in the output sparse, symmetric matrix `M` is given by:
+// Each entry in the output sparse, symmetric matrix `M` is given by:
 //
 // \\[
 // M_{ij} = \begin{cases}
diff --git a/include/smooth.h b/include/smooth.h
index 6765fe2..38fad4b 100644
--- a/include/smooth.h
+++ b/include/smooth.h
@@ -2,6 +2,9 @@
 #define SMOOTH_H
 #include <Eigen/Core>
 #include <Eigen/Sparse>
+#include <igl/edge_lengths.h>
+#include "cotmatrix.h"
+#include "massmatrix.h"
 // Given a mesh (`V`,`F`) and data specified per-vertex (`G`), smooth this data
 // using a single implicit Laplacian smoothing step.
 //
diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..f3d4e54 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,10 +1,49 @@
 #include "cotmatrix.h"
 
+// Calcutes 1/2 * cot(theta) in a triangle with side lengths
+// a, b, and c, where theta is the angle opposite c.
+double halfcot(
+    double a,
+    double b, 
+    double c) 
+{
+    // Heron's formula for area
+    double s = (a + b + c) / 2;
+    double A = std::sqrt(s * (s - a) * (s - b) * (s - c));
+
+    // cos(theta) / sin(theta)
+    return (std::pow(c, 2) - std::pow(a, 2) - std::pow(b, 2)) / (-8 * A);
+}
+
 void cotmatrix(
-  const Eigen::MatrixXd & l,
-  const Eigen::MatrixXi & F,
-  Eigen::SparseMatrix<double> & L)
+    const Eigen::MatrixXd & l,
+    const Eigen::MatrixXi & F,
+    Eigen::SparseMatrix<double> & L)
 {
-  // Add your code here
+    typedef Eigen::Triplet<double> T;
+    std::vector<T> triplets;
+    triplets.reserve(F.size() * 4);
+
+    for (int i = 0; i < F.rows(); i++) 
+    {
+        Eigen::Vector3i f = F.row(i);
+        Eigen::Vector3d lengths = l.row(i);
+
+        for (int j = 0; j < f.size(); j++) 
+        {
+            // lengths[j] corresponds to the length of the edge 
+            // between vertices f[(j + 1) % 3] and f[(j + 2) % 3]
+            double c = halfcot(lengths[(j + 1) % 3], lengths[(j + 2) % 3], lengths[j]);
+            triplets.push_back(T(f[(j + 1) % 3], f[(j + 2) % 3], c));
+            triplets.push_back(T(f[(j + 2) % 3], f[(j + 1) % 3], c));
+
+            // populate diagonal with negative sum of off-diagonal entries
+            triplets.push_back(T(f[(j + 1) % 3], f[(j + 1) % 3], -c));
+            triplets.push_back(T(f[(j + 2) % 3], f[(j + 2) % 3], -c));
+        }
+    }
+
+    L.resize(F.maxCoeff() + 1, F.maxCoeff() + 1);
+    L.setFromTriplets(triplets.begin(), triplets.end());
 }
 
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..514568d 100644
--- a/src/massmatrix.cpp
+++ b/src/massmatrix.cpp
@@ -1,10 +1,29 @@
 #include "massmatrix.h"
 
 void massmatrix(
-  const Eigen::MatrixXd & l,
-  const Eigen::MatrixXi & F,
-  Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
+    const Eigen::MatrixXd & l,
+    const Eigen::MatrixXi & F,
+    Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
 {
-  // Add your code here
+    Eigen::VectorXd dblA(F.rows());
+    igl::doublearea(l, 0, dblA);
+
+    Eigen::VectorXd areas = Eigen::VectorXd::Zero(F.maxCoeff() + 1);
+    
+    for (int i = 0; i < F.rows(); i++) 
+    {
+        Eigen::Vector3i f = F.row(i);
+        double area = dblA[i];
+
+        for (int j = 0; j < f.size(); j++) 
+        {
+            int v = f[j];
+            areas[v] += area;
+
+        }
+    }
+
+    M.resize(F.maxCoeff() + 1);
+    M.diagonal() = 1 / (double)6 * areas;
 }
 
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..e93b19f 100644
--- a/src/smooth.cpp
+++ b/src/smooth.cpp
@@ -7,6 +7,24 @@ void smooth(
     double lambda,
     Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = G;
+    Eigen::MatrixXd l;
+    igl::edge_lengths(V, F, l);
+
+    Eigen::DiagonalMatrix<double, Eigen::Dynamic> M;
+    massmatrix(l, F, M);
+
+    Eigen::SparseMatrix<double> L;
+    cotmatrix(l, F, L);
+
+    Eigen::MatrixXd MG = M * G;
+
+    // A = M - lamda * L
+    Eigen::SparseMatrix<double> A = -lambda * L;
+    for (int i = 0; i < A.rows(); i++) 
+    {
+        A.coeffRef(i, i) += M.diagonal()[i];
+    }
+
+    Eigen::SimplicialCholesky<Eigen::SparseMatrix<double>> chol(A);
+    U = chol.solve(MG);
 }