diff --git a/src/lscm.cpp b/src/lscm.cpp
index d6e1b6c..45209de 100644
--- a/src/lscm.cpp
+++ b/src/lscm.cpp
@@ -1,10 +1,40 @@
 #include "lscm.h"
+#include <igl/eigs.h>
+#include <igl/repdiag.h>
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include "vector_area_matrix.h"
+#include <Eigen/SVD>
+
+using namespace Eigen;
+using namespace std;
 
 void lscm(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = V.leftCols(2);
+  // construct matrix Q
+  SparseMatrix<double> L, L_rep, Q, A, M, B;
+  vector_area_matrix(F, A);
+  igl::cotmatrix(V, F, L);
+  igl::repdiag(L, 2, L_rep);
+  Q = L_rep - A;
+
+  // similarly, construct matrix B
+  igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_DEFAULT, M);
+  igl::repdiag(M, 2, B);
+
+  // start to solve the generalized Eigen value problem
+  MatrixXd eVector;
+  VectorXd eValue;
+  igl::eigs(Q, B, 3, igl::EIGS_TYPE_SM, eVector, eValue);
+  // extract the Fiedler vector: 3rd -- first 2 trivial
+  int n = V.rows();
+  U.resize(n, 2);
+  U.col(0) = eVector.col(2).block(0, 0, n, 1);
+  U.col(1) = eVector.col(2).block(n, 0, n, 1);
+  // using PCA: SVD solves
+  JacobiSVD<MatrixXd> SVD(U.transpose() * U, ComputeThinU | ComputeThinV);
+  U = U * SVD.matrixU();
 }
diff --git a/src/tutte.cpp b/src/tutte.cpp
index 6f6109b..957cbd9 100644
--- a/src/tutte.cpp
+++ b/src/tutte.cpp
@@ -1,11 +1,31 @@
 #include "tutte.h"
+#include <igl/cotmatrix.h>
+#include <igl/boundary_loop.h>
+#include <igl/map_vertices_to_circle.h>
+#include <igl/min_quad_with_fixed.h>
+
+using namespace Eigen;
 
 void tutte(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = V.leftCols(2);
-}
+  // compute the cotangent matrix L
+  SparseMatrix<double> L;
+  igl::cotmatrix(V, F, L);
 
+  // construct boundary constraints: convex polygon UV
+  VectorXi b;
+  igl::boundary_loop(F, b);
+  MatrixXd UV;
+  igl::map_vertices_to_circle(V, b, UV);
+
+  // compute U
+  SparseMatrix<double> Aeq;
+  igl::min_quad_with_fixed_data<double> data;
+  igl::min_quad_with_fixed_precompute(L, b, Aeq, false, data);
+  MatrixXd Beq, B;
+  B = MatrixXd::Zero(V.rows(), 2);
+  igl::min_quad_with_fixed_solve(data, B, UV, Beq, U);
+}
diff --git a/src/vector_area_matrix.cpp b/src/vector_area_matrix.cpp
index 1d60ef8..f68e39f 100644
--- a/src/vector_area_matrix.cpp
+++ b/src/vector_area_matrix.cpp
@@ -1,11 +1,42 @@
 #include "vector_area_matrix.h"
+#include <igl/boundary_loop.h>
+
+using namespace Eigen;
+using namespace std;
 
 void vector_area_matrix(
   const Eigen::MatrixXi & F,
   Eigen::SparseMatrix<double>& A)
 {
-  // Replace with your code
-  int V_size = F.maxCoeff()+1;
-  A.resize(V_size*2,V_size*2);
-}
+  int n = F.maxCoeff() + 1;
+  A.resize(n * 2, n * 2);
+
+  vector<vector<int>> b;
+  igl::boundary_loop(F, b);
 
+  // construct A
+  // two points scenario
+  //  0   0   0  0.5
+  //  0   0 -0.5   0
+  //  0 -0.5   0   0
+  // 0.5   0   0   0
+  vector<Triplet<double>> triplets;
+  for (int i = 0; i < b.size(); i++) {
+    vector<int> bd = b[i];
+    for (int j = 0; j< bd.size(); j++) {
+      int p1 = bd[j];
+      int p2_index = (j + 1) % bd.size();
+      int p2 = bd[p2_index];
+
+      // ensure symmetry
+      // A
+      triplets.push_back(Triplet<double>(p1, p2 + n, 0.5));
+      triplets.push_back(Triplet<double>(p2, p1 + n, -0.5));
+      // A_T
+      triplets.push_back(Triplet<double>(p2 + n, p1, 0.5));
+      triplets.push_back(Triplet<double>(p1 + n, p2, -0.5));
+    }
+  }
+  A.setFromTriplets(triplets.begin(), triplets.end());
+  A = 0.5 * A;
+}