diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..41937d0 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,10 +1,37 @@
 #include "cotmatrix.h"
+#include <cmath>
+
+typedef Eigen::Triplet<double> T;
 
 void cotmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::SparseMatrix<double> & L)
 {
-  // Add your code here
+    std::vector<Eigen::Triplet<double>> trips;
+    L.resize(F.maxCoeff()+1, F.maxCoeff()+1);
+    
+    for (int i = 0; i < F.rows(); i++) {
+        for (int j_ = 0; j_ < 3; j_++) {
+            int j[] = {j_, (j_+1)%3, (j_+2)%3};
+            double a = l(i,j[2]);
+            double b = l(i,j[0]);
+            double c = l(i,j[1]);
+            
+            //Heron's formula to calculate area
+            double s = (a+b+c) / 2.0;
+            double A = sqrt(s * (s-a) * (s-b) * (s-c));
+            //half the cotangent = cos / sin / 2
+            double cot = (((a*a-b*b-c*c)/(-2*b*c)) / ((2*A)/(b*c))) / 2;
+            
+            //put to matrix
+            trips.push_back(T(F(i,j[0]), F(i,j[1]), cot));
+            trips.push_back(T(F(i,j[1]), F(i,j[0]), cot));
+            trips.push_back(T(F(i,j[0]), F(i,j[0]), -cot));
+            trips.push_back(T(F(i,j[1]), F(i,j[1]), -cot));
+        }
+    }
+    
+    L.setFromTriplets(trips.begin(), trips.end());
 }
 
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..f7c0af5 100644
--- a/src/massmatrix.cpp
+++ b/src/massmatrix.cpp
@@ -1,10 +1,20 @@
 #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);
+    M.setZero();
+    Eigen::MatrixXd area;
+    //compute double the area
+    igl::doublearea(l, area);
+    for (int i = 0; i < F.rows(); i++) {
+        //add 1/3 of the area to each
+        M.diagonal()(F(i, 0)) += area(i) / 2.0 / 3;
+        M.diagonal()(F(i, 1)) += area(i) / 2.0 / 3;
+        M.diagonal()(F(i, 2)) += area(i) / 2.0 / 3;
+    }
 }
-
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..b046e36 100644
--- a/src/smooth.cpp
+++ b/src/smooth.cpp
@@ -1,4 +1,8 @@
 #include "smooth.h"
+#include "massmatrix.h"
+#include "cotmatrix.h"
+#include <Eigen/SparseCholesky>
+#include <igl/edge_lengths.h>
 
 void smooth(
     const Eigen::MatrixXd & V,
@@ -7,6 +11,22 @@ void smooth(
     double lambda,
     Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = G;
+    //edge lengths, mass matrix, cotangent matrix
+    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::SimplicialLDLT<Eigen::SparseMatrix<double>> S;
+    
+    //A = M - (lambda)*L
+    Eigen::SparseMatrix<double> A = -lambda * L;
+    for (int i = 0; i < M.rows(); i++) {
+        A.coeffRef(i, i) += M.diagonal()[i];
+    }
+    
+    S.compute(A);
+    //LDLT solve
+    U = S.solve(M * G);
 }