diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..56fae68 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,4 +1,6 @@
 #include "cotmatrix.h"
+#include <iostream>
+#include <math.h>
 
 void cotmatrix(
   const Eigen::MatrixXd & l,
@@ -6,5 +8,45 @@ void cotmatrix(
   Eigen::SparseMatrix<double> & L)
 {
   // Add your code here
+	L.resize(F.maxCoeff() + 1, F.maxCoeff() + 1);
+	L.setZero();
+	std::vector<Eigen::Triplet<double>> triplets;
+	triplets.reserve(6 * F.rows());
+
+	for (int f = 0; f < F.rows(); f++) {
+		int i = F(f, 0); int j = F(f, 1); int k = F(f, 2);
+		double a = l(f, 0); double b = l(f, 1); double c = l(f, 2);
+		double s = (a + b + c) / 2.0;
+		double A = sqrt(s * (s-a) * (s-b) * (s-c));
+		double cot;
+		// Consider edge ij
+		cot = (pow(a, 2) + pow(b, 2) - pow(c, 2)) / (4 * A);
+		triplets.push_back({ i, j, 0.5 * cot });
+		triplets.push_back({ j, i, 0.5 * cot });
+
+		// Consider edge jk
+		cot = (pow(b, 2) + pow(c, 2) - pow(a, 2)) / (4 * A);
+		triplets.push_back({ j, k, 0.5 * cot });
+		triplets.push_back({ k, j, 0.5 * cot });
+
+		// Consider edge ik
+		cot = (pow(a, 2) + pow(c, 2) - pow(b, 2)) / (4 * A);
+		triplets.push_back({ i, k, 0.5 * cot });
+		triplets.push_back({ k, i, 0.5 * cot });
+	}
+	L.setFromTriplets(triplets.begin(), triplets.end());
+
+	// Fill in the diagonal entries
+	Eigen::SparseMatrix<double> Ld(F.maxCoeff() + 1, F.maxCoeff() + 1);
+	Ld.setZero();
+	std::vector<Eigen::Triplet<double>> tripletsd;
+	tripletsd.reserve(L.rows());
+	for (int f = 0; f < L.rows(); f++) {
+		tripletsd.push_back({ f, f,  L.row(f).sum() });
+	}
+
+	Ld.setFromTriplets(tripletsd.begin(), tripletsd.end());
+
+	L = L - Ld;
 }
 
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..51fdd55 100644
--- a/src/massmatrix.cpp
+++ b/src/massmatrix.cpp
@@ -1,4 +1,6 @@
 #include "massmatrix.h"
+#include <iostream>
+#include <igl/doublearea.h>
 
 void massmatrix(
   const Eigen::MatrixXd & l,
@@ -6,5 +8,16 @@ void massmatrix(
   Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
 {
   // Add your code here
-}
+	int nv = F.maxCoeff() + 1;
+	M.resize(F.maxCoeff() + 1);
+	M.setZero();
+	Eigen::MatrixXd dblA(F.rows(), 1);
+	igl::doublearea(l, dblA);
 
+	// Fill in the diagonal entries of M
+	for (int f = 0; f < F.rows(); f++) {
+		M.diagonal()[F(f, 0)] += (1.0 / 6.0) * dblA(f);
+		M.diagonal()[F(f, 1)] += (1.0 / 6.0) * dblA(f);
+		M.diagonal()[F(f, 2)] += (1.0 / 6.0) * dblA(f);
+	}
+}
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..621f28f 100644
--- a/src/smooth.cpp
+++ b/src/smooth.cpp
@@ -1,4 +1,11 @@
 #include "smooth.h"
+#include <iostream>
+#include <igl/edge_lengths.h>
+#include <igl/doublearea.h>
+#include "cotmatrix.h"
+#include "massmatrix.h"
+#include<Eigen/SparseCholesky>
+#include<Eigen/IterativeLinearSolvers>
 
 void smooth(
     const Eigen::MatrixXd & V,
@@ -8,5 +15,37 @@ void smooth(
     Eigen::MatrixXd & U)
 {
   // Replace with your code
-  U = G;
+
+  // Calculate edge lengths
+  Eigen::MatrixXd l;
+  igl::edge_lengths(V, F, l);
+
+  // Calculate cotmatrix
+  Eigen::SparseMatrix<double> L;
+  cotmatrix(l, F, L);
+  
+  // Calculate massmatrix
+  Eigen::DiagonalMatrix<double, Eigen::Dynamic> M;
+  massmatrix(l, F, M);
+
+  // Calculate M - lambda * L
+  Eigen::SparseMatrix<double> A;
+  Eigen::SparseMatrix<double> Md(F.maxCoeff() + 1, F.maxCoeff() + 1);
+  Md.setZero();
+  std::vector<Eigen::Triplet<double>> tripletsd;
+  tripletsd.reserve(F.rows());
+  for (int f = 0; f < L.rows(); f++) {
+	  tripletsd.push_back({ f, f,  M.diagonal()[f] });
+  }
+  Md.setFromTriplets(tripletsd.begin(), tripletsd.end());
+  A.resize(L.rows(), L.cols());
+  A = Md - lambda * L;
+
+  // Solve the linear system
+  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> sparseSolver;
+  sparseSolver.compute(A);
+  U = sparseSolver.solve(M * G);
+
+
+
 }