diff --git a/src/lscm.cpp b/src/lscm.cpp
index d6e1b6c..b8eb29c 100644
--- a/src/lscm.cpp
+++ b/src/lscm.cpp
@@ -1,10 +1,41 @@
 #include "lscm.h"
+#include <Eigen/SparseCore>
+#include "vector_area_matrix.h"
+#include <igl/cotmatrix.h>
+#include <igl/repdiag.h>
+#include <igl/massmatrix.h>
+#include <igl/eigs.h>
+#include <Eigen/SVD>
+#include <Eigen/LU>
 
 void lscm(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  Eigen::MatrixXd & U)
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    Eigen::MatrixXd & U)
 {
   // Replace with your code
-  U = V.leftCols(2);
+  int n = V.rows();
+  U.resize(n, 2);
+
+  Eigen::SparseMatrix<double> A, Q, L, L2;
+
+  vector_area_matrix(F, A);
+  igl::cotmatrix(V, F, L);
+  igl::repdiag(L, 2, L2);
+  Q = L2 - A;
+
+  Eigen::SparseMatrix<double> M, B;
+  igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_DEFAULT, M);
+  igl::repdiag(M, 2, B);
+
+  // Solve general eigen value problem
+  Eigen::MatrixXd eigenVec;
+  Eigen::VectorXd eigenVal;
+  igl::eigs(Q, B, 3, igl::EIGS_TYPE_SM, eigenVec, eigenVal);
+  U.col(0) << eigenVec.col(2).head(n);
+  U.col(1) << eigenVec.col(2).tail(n);
+
+  // SVD
+  Eigen::JacobiSVD<Eigen::MatrixXd> svd(U.transpose() * U, Eigen::ComputeFullU | Eigen::ComputeFullV);
+  U = U * svd.matrixU();
 }
diff --git a/src/tutte.cpp b/src/tutte.cpp
index 6f6109b..291ab78 100644
--- a/src/tutte.cpp
+++ b/src/tutte.cpp
@@ -1,11 +1,29 @@
 #include "tutte.h"
+#include <Eigen/Sparse>
+#include <igl/boundary_loop.h>
+#include <igl/cotmatrix.h>
+#include <igl/map_vertices_to_circle.h>
+#include <igl/min_quad_with_fixed.h>
 
 void tutte(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = V.leftCols(2);
+  Eigen::VectorXi bL;
+  igl::boundary_loop(F, bL);
+
+  Eigen::MatrixXd UV;
+  igl::map_vertices_to_circle(V, bL, UV);
+
+  Eigen::SparseMatrix<double> L(V.rows(), V.rows());
+  igl::cotmatrix(V, F, L);
+
+  igl::min_quad_with_fixed_data<double> data;
+  igl::min_quad_with_fixed_precompute(L, bL, Eigen::SparseMatrix<double>(), false, data);
+
+  Eigen::VectorXd B = Eigen::VectorXd::Zero(data.n, 1);
+  igl::min_quad_with_fixed_solve(data, B, UV, Eigen::MatrixXd(), U);
+  U.col(0) = -U.col(0);
 }
 
diff --git a/src/vector_area_matrix.cpp b/src/vector_area_matrix.cpp
index 1d60ef8..775a740 100644
--- a/src/vector_area_matrix.cpp
+++ b/src/vector_area_matrix.cpp
@@ -1,11 +1,29 @@
 #include "vector_area_matrix.h"
+#include <igl/boundary_loop.h>
 
 void vector_area_matrix(
-  const Eigen::MatrixXi & F,
-  Eigen::SparseMatrix<double>& A)
+    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);
+
+  std::vector<std::vector<int>> bL;
+  igl::boundary_loop(F, bL);
+
+  for (int x = 0; x < bL.size(); x++) {
+    std::vector<int> b = bL[x];
+    for (int y = 0; y < b.size(); y++) {
+      int i = b[y];
+      int j = b[(y + 1) % b.size()];
+      A.coeffRef(i, n + j) = 0.5;
+      A.coeffRef(j, n + i) = -0.5;
+    }
+  }
+
+  // A to be symmetric
+  Eigen::SparseMatrix<double> AT = A.transpose();
+  A = A + AT;
+  // However A *= 0.5 gives strange pattern.
 }