diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..720cc30 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,4 +1,5 @@
 #include "cotmatrix.h"
+#include <igl/doublearea.h>
 
 void cotmatrix(
   const Eigen::MatrixXd & l,
@@ -6,5 +7,39 @@ void cotmatrix(
   Eigen::SparseMatrix<double> & L)
 {
   // Add your code here
+
+	// compute area for all facets
+	Eigen::VectorXd doublearea(F.rows());
+	igl::doublearea(l, doublearea);
+
+	// compute 1/2*cot for angle inside each facet
+	typedef Eigen::Triplet<double> T;
+	typedef Eigen::SparseMatrix<double> SpMat;
+	std::vector<T> triplets;
+	triplets.reserve(F.rows() * 4);
+	for (int i = 0; i < F.rows(); i++) {
+		double a_2 = l(i, 2)*l(i, 2);
+		double b_2 = l(i, 0)*l(i, 0);
+		double c_2 = l(i, 1)*l(i, 1);
+		double douArea_inv = 1 / doublearea(i);
+		double cot_a = 0.5*(b_2 + c_2 - a_2) * douArea_inv;
+		double cot_b = 0.5*(a_2 + c_2 - b_2) * douArea_inv;
+		double cot_c = 0.5*(a_2 + b_2 - c_2) * douArea_inv;
+		triplets.push_back(T(F(i, 0), F(i, 1), 0.5*cot_a));
+		triplets.push_back(T(F(i, 1), F(i, 2), 0.5*cot_b));
+		triplets.push_back(T(F(i, 2), F(i, 0), 0.5*cot_c));
+	}
+	int V_count = F.maxCoeff() + 1;
+	SpMat cotMat(V_count, V_count);
+	cotMat.setFromTriplets(triplets.begin(), triplets.end());
+
+	// compute 1/2*cot alpha_ij + 1/2*cot beta_ij for each edge ij that exists
+	SpMat cotMat_t = cotMat.transpose();
+	L = cotMat_t + cotMat;
+
+	// compute diagonal part
+	for (int i = 0; i < V_count; i++) {
+		L.insert(i, i) = -L.row(i).sum();
+	}
 }
 
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..43c2be3 100644
--- a/src/massmatrix.cpp
+++ b/src/massmatrix.cpp
@@ -1,4 +1,5 @@
 #include "massmatrix.h"
+#include <igl/doublearea.h>
 
 void massmatrix(
   const Eigen::MatrixXd & l,
@@ -6,5 +7,20 @@ void massmatrix(
   Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
 {
   // Add your code here
+	
+	// compute area for all facets
+	Eigen::VectorXd doublearea(F.rows());
+	igl::doublearea(l, doublearea);
+	Eigen::VectorXd area = 0.5*doublearea;
+	// compute diagonal of M as a vector
+	int V_count = F.maxCoeff() + 1;
+	Eigen::VectorXd diagonal = Eigen::VectorXd::Zero(V_count);
+	for (int i = 0; i < F.rows(); i++) {
+		diagonal(F(i, 0)) += area(i);
+		diagonal(F(i, 1)) += area(i);
+		diagonal(F(i, 2)) += area(i);
+	}
+	diagonal = diagonal / 3.0;
+	M.diagonal() = diagonal;
 }
 
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..381781c 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,
@@ -8,5 +11,27 @@ void smooth(
     Eigen::MatrixXd & U)
 {
   // Replace with your code
-  U = G;
+	// compute edge length
+	Eigen::MatrixXd l;
+	igl::edge_lengths(V, F, l);
+
+	// compute L and M
+	Eigen::SparseMatrix<double> L;
+	Eigen::DiagonalMatrix<double, Eigen::Dynamic> M;
+	cotmatrix(l, F, L);
+	massmatrix(l, F, M);
+
+	// solve U
+	Eigen::VectorXd diagM = M.diagonal();
+	typedef Eigen::Triplet<double> T;
+	typedef Eigen::SparseMatrix<double> SpMat;
+	std::vector<T> triplets;
+	triplets.reserve(diagM.size()*2);
+	for (int i = 0; i < diagM.size(); i++) {
+		triplets.push_back(T(i, i, diagM(i)));
+	}
+	SpMat spaM(M.rows(), M.cols());
+	spaM.setFromTriplets(triplets.begin(), triplets.end());
+	Eigen::SimplicialCholesky<Eigen::SparseMatrix<double>> chol(spaM - lambda * L);
+	U = chol.solve(M*G);
 }