diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..059ee64 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,10 +1,55 @@
 #include "cotmatrix.h"
+#include <iostream>
+#include <vector>
+#include <cmath>
+
+typedef Eigen::Triplet<double> T;
+
+// Given edge lengths of a triangle ABC
+// Return the cotan of angle A
+double my_cot(double a, double b, double c) {
+	// Heron's Formula
+	double s = (a + b + c) / 2.0;
+	double area = sqrt(s * (s - a) * (s - b) * (s - c));
+
+	// Law of Sine and Cosine
+	double sinA = 2.0 * area / (b * c);
+	double cosA = (a*a - b*b - c*c) / (-2.0 * b * c);
+
+	return cosA / sinA;
+}
 
 void cotmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::SparseMatrix<double> & L)
 {
-  // Add your code here
+	int nV = F.maxCoeff() + 1;
+	std::vector<T> tripletList;
+
+	for(int fi = 0; fi < F.rows(); fi++) {
+		for(int e = 0; e < 3; e++) {
+			int i = F(fi, e); 
+			int j = F(fi, (e + 1) % 3);
+
+			double a = l(fi, (e + 2) % 3);
+			double b = l(fi, (e + 0) % 3);
+			double c = l(fi, (e + 1) % 3);
+
+			double cotA = my_cot(a, b, c);
+
+			double Lij = 0.5 * cotA;
+
+			tripletList.push_back(T(i, j, Lij)); // Duplicates will be summed
+			tripletList.push_back(T(j, i, Lij));
+
+			// Diagonal
+			tripletList.push_back(T(i, i, -Lij));
+			tripletList.push_back(T(j, j, -Lij));
+		}
+	}
+
+	L.resize(nV, nV);
+	L.setFromTriplets(tripletList.begin(), tripletList.end());
 }
 
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..1a6b3b3 100644
--- a/src/massmatrix.cpp
+++ b/src/massmatrix.cpp
@@ -1,10 +1,22 @@
 #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
+  int nV = F.maxCoeff() + 1;
+  M.resize(nV);
+
+  Eigen::VectorXd dblA;
+  igl::doublearea(l, dblA);
+
+  for(int fi = 0; fi < F.rows(); fi++) {
+  	for(int vd = 0; vd < 3; vd++) {
+  		int vert = F(fi, vd);
+  		M.diagonal()[vert] += dblA[fi] / 6.0;
+  	}
+  }
 }
 
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..87a7c57 100644
--- a/src/smooth.cpp
+++ b/src/smooth.cpp
@@ -1,4 +1,7 @@
 #include "smooth.h"
+#include "cotmatrix.h"
+#include "massmatrix.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);
+
+  // M * G = (M - lambda * L) * U 
+  Eigen::SparseMatrix<double> A = -lambda * L;
+  for(int i = 0; i < M.rows(); i++) {
+  	A.coeffRef(i, i) += M.diagonal()[i];
+  }
+
+  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > solver(A);
+  U = solver.solve(M * G);
 }