diff --git a/src/lscm.cpp b/src/lscm.cpp
index d6e1b6c..42550af 100644
--- a/src/lscm.cpp
+++ b/src/lscm.cpp
@@ -1,10 +1,36 @@
 #include "lscm.h"
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <igl/eigs.h>
+#include <igl/repdiag.h>
+#include "vector_area_matrix.h"
 
 void lscm(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = V.leftCols(2);
+  U.resize(V.rows(), 2);
+  // Quardratic Coefficients:
+  Eigen::SparseMatrix<double> L, Q, A;
+  igl::cotmatrix(V, F, L);
+  vector_area_matrix(F, A);
+  igl::repdiag(L, 2, Q);
+  Q = Q - A;
+
+  // Free Boundary:
+  Eigen::SparseMatrix<double> M, B;
+  igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_DEFAULT, M);
+  igl::repdiag(M, 2, B);
+  Eigen::MatrixXd sU;
+  Eigen::VectorXd sS;
+  igl::eigs(Q, B, 3, igl::EIGS_TYPE_SM, sU, sS);
+  // First two are trivial solutions
+  U.col(0) = sU.col(2).head(V.rows());
+  U.col(1) = sU.col(2).tail(V.rows());
+
+  // Apply 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..182c284 100644
--- a/src/tutte.cpp
+++ b/src/tutte.cpp
@@ -1,11 +1,32 @@
 #include "tutte.h"
+#include <igl/boundary_loop.h>
+#include <igl/map_vertices_to_circle.h>
+#include <igl/min_quad_with_fixed.h>
+#include <igl/cotmatrix.h>
 
 void tutte(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = V.leftCols(2);
+  // Known coefficients: Boundary indexes and mapped vertices
+  Eigen::VectorXi bound_index;
+  Eigen::MatrixXd bound_vertices;
+  igl::boundary_loop(F, bound_index);
+  igl::map_vertices_to_circle(V, bound_index, bound_vertices);
+
+  // Quadratic coefficients matrix: Contagent Laplacian
+  Eigen::SparseMatrix<double> L;
+  igl::cotmatrix(V, F, L);
+
+  // Precompute quadratic functions:
+  igl::min_quad_with_fixed_data<double> data;
+  igl::min_quad_with_fixed_precompute(L, bound_index, 
+                Eigen::SparseMatrix<double>(), false, data);
+  
+  // Solve quardratic functions:
+  Eigen::MatrixXd B = Eigen::MatrixXd::Zero(data.n, 2);
+  igl::min_quad_with_fixed_solve(data, B, bound_vertices, 
+                Eigen::MatrixXd(), U);
 }
 
diff --git a/src/vector_area_matrix.cpp b/src/vector_area_matrix.cpp
index 1d60ef8..01740fc 100644
--- a/src/vector_area_matrix.cpp
+++ b/src/vector_area_matrix.cpp
@@ -1,11 +1,42 @@
 #include "vector_area_matrix.h"
+#include <Eigen/Eigen>
+#include <igl/boundary_loop.h>
+#include <igl/vector_area_matrix.h>
+typedef Eigen::Triplet<double> T;
 
 void vector_area_matrix(
   const Eigen::MatrixXi & F,
   Eigen::SparseMatrix<double>& A)
 {
-  // Replace with your code
   int V_size = F.maxCoeff()+1;
+  int n = V_size;
   A.resize(V_size*2,V_size*2);
+
+  // Get boundary loop vertices:
+  std::vector<std::vector<int>> bound_loop;
+  igl::boundary_loop(F, bound_loop);
+
+  // Calculate area matrix:
+  int loop_num = bound_loop.size();
+  std::vector<T> triplet_list;
+
+  for (int i = 0; i < loop_num; i ++) {
+    std::vector<int> cur_loop = bound_loop[i];
+    int n = cur_loop.size();
+    for (int j = 0; j < cur_loop.size(); j ++) {
+      int ui = cur_loop[j];
+      int uj = cur_loop[(j + 1) % n];
+
+      // det(ui, uj):
+      // uiu * ujv - uiv * uju
+      triplet_list.push_back(T(ui, uj + V_size, 0.5));
+      triplet_list.push_back(T(ui + V_size, uj, -0.5));
+    }
+  }
+
+  // Symmetric:
+  A.setFromTriplets(triplet_list.begin(), triplet_list.end());
+  Eigen::SparseMatrix<double> AT =  A.transpose();
+  A = A + AT; // times 0.5 here seems incorrect
 }