diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..ea3d1c9 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,4 +1,20 @@
 #include "cotmatrix.h"
+#include <cmath>
+#include <vector>
+#include <iostream>
+//calculate cos C
+double cos_abc(double a, double b, double c){
+    return (a * a + b * b - c * c) / (2 * a * b);
+}
+
+//calculate sin C
+double sin_abc(double a, double b, double c){
+    double s = (a + b + c)/2;
+    double area = sqrt(s * (s - a) * (s -b) * (s - c));
+    //return c / (a*b*c/(2 * area));
+    return 2 * area/(a * b);
+}
+
 
 void cotmatrix(
   const Eigen::MatrixXd & l,
@@ -6,5 +22,33 @@ void cotmatrix(
   Eigen::SparseMatrix<double> & L)
 {
   // Add your code here
+    L.resize(F.maxCoeff() + 1, F.maxCoeff() + 1);
+    for (int row = 0; row < F.rows(); ++row){
+        for (int col = 0; col < F.cols(); ++col){ //3 edges [0, 1], [1, 2] [2, 0]
+            int i = F(row, col);
+            int j = F(row, (col + 1) % F.cols());
+            double a = l(row, col);
+            double b = l(row, (col +1) % l.cols());
+            double c = l(row, (col +2) % l.cols());
+            L.coeffRef(i,j) += 0.5 * cos_abc(a, b, c)/sin_abc(a, b, c);
+            //std::cout <<i << " " << j << " " <<0.5 * cos_abc(a, b, c)/sin_abc(a, b, c) << std::endl;
+        } 
+        for (int col = 0; col < F.cols(); ++col){ //3 edges [1, 0], [2, 1] [0, 2] symmetric
+            int j = F(row, col);
+            int i = F(row, (col + 1) % F.cols());
+            L.coeffRef(i,j) = L.coeffRef(j,i);
+            //std::cout <<i << " " << j << " " <<0.5 * cos_abc(a, b, c)/sin_abc(a, b, c) << std::endl;
+        }
+        for (int col = 0; col < F.cols(); ++col){ //diag
+            int i = F(row, col);
+            double a = l(row, col); //[1,2]
+            double b = l(row, (col +1) % l.cols()); //[0,2]
+            double c = l(row, (col +2) % l.cols()); //[0,1]
+            L.coeffRef(i,i) -= 0.5 * cos_abc(a, b, c)/sin_abc(a, b, c);
+            
+            L.coeffRef(i,i) -= 0.5 * cos_abc(a, c, b)/sin_abc(a, c, b);
+
+        }
+    }
 }
 
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..1ac0aac 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
+    Eigen::MatrixXd areas;
+    igl::doublearea(l, areas);
+    M.resize(F.maxCoeff() + 1);
+    for (int row = 0; row < F.rows(); ++row){
+        for (int col = 0; col < F.cols(); ++col){
+            int i = F(row, col);
+            M.diagonal()[i] += 1/3.0 * 0.5 * areas(row);
+        }
+    }
+
 }
 
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..cdde68d 100644
--- a/src/smooth.cpp
+++ b/src/smooth.cpp
@@ -1,5 +1,8 @@
 #include "smooth.h"
-
+#include "cotmatrix.h"
+#include "massmatrix.h"
+#include <igl/edge_lengths.h>
+#include <iostream>
 void smooth(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
@@ -8,5 +11,22 @@ void smooth(
     Eigen::MatrixXd & U)
 {
   // Replace with your code
-  U = G;
+  Eigen::DiagonalMatrix<double,Eigen::Dynamic>  M;
+  Eigen::MatrixXd l;
+  Eigen::SparseMatrix<double> L;
+  //std::cout << "1" << std::endl;
+  igl::edge_lengths(V, F, l);
+  cotmatrix(l, F, L);
+  //std::cout << "2" << std::endl;
+  massmatrix(l, F, M);
+  //std::cout << "3" << std::endl;
+  Eigen::MatrixXd left = M * G;
+  Eigen::MatrixXd temp(M);
+  Eigen::SparseMatrix<double> right =  temp.sparseView();
+  right = right - lambda * L;
+  //std::cout << "4" << std::endl;
+  Eigen::BiCGSTAB<Eigen::SparseMatrix<double>> cg;
+  cg.compute(right);
+  U = cg.solve(left);
+
 }