diff --git a/src/lscm.cpp b/src/lscm.cpp
index d6e1b6c..ebb41ed 100644
--- a/src/lscm.cpp
+++ b/src/lscm.cpp
@@ -1,4 +1,15 @@
 #include "lscm.h"
+#include "igl/cotmatrix.h"
+#include "igl/massmatrix.h"
+#include "igl/repdiag.h"
+#include "igl/eigs.h"
+#include "Eigen/src/SVD/JacobiSVD.h"
+#include "vector_area_matrix.h"
+#include <iostream>
+
+using namespace std;
+using namespace igl;
+using namespace Eigen;
 
 void lscm(
   const Eigen::MatrixXd & V,
@@ -7,4 +18,30 @@ void lscm(
 {
   // Replace with your code
   U = V.leftCols(2);
+  // Calculate L
+  Eigen::SparseMatrix<double> L(V.rows(), V.rows());
+  L.setZero();
+  igl::cotmatrix(V, F, L);
+  Eigen::SparseMatrix<double> A;
+  vector_area_matrix(F, A);
+
+  // Set up the quadratic matrices
+  Eigen::SparseMatrix<double> L_rep;
+  igl::repdiag(L, 2, L_rep);
+  Eigen::SparseMatrix<double> Q = L_rep - A;
+  Eigen::SparseMatrix<double> M;
+  igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_DEFAULT, M);
+  Eigen::SparseMatrix<double> B;
+  igl::repdiag(M, 2, B);
+
+  // Solving the generalized eigen vector problem
+  Eigen::MatrixXd sU;
+  Eigen::VectorXd sS;
+  igl::eigs(Q, B, 3, igl::EIGS_TYPE_SM, sU, sS);
+  U.col(0) = sU.block(0, 2, V.rows(), 1);
+  U.col(1) = sU.block(V.rows(), 2, V.rows(), 1);
+
+  // Align the parameterization
+  Eigen::JacobiSVD<Eigen::MatrixXd> svd(U.transpose() * U, Eigen::ComputeThinU | Eigen::ComputeThinV);
+  U = U * svd.matrixV();
 }
diff --git a/src/tutte.cpp b/src/tutte.cpp
index 6f6109b..8cf6ba1 100644
--- a/src/tutte.cpp
+++ b/src/tutte.cpp
@@ -1,4 +1,11 @@
 #include "tutte.h"
+#include <iostream>
+#include <igl/map_vertices_to_circle.h>
+#include <igl/boundary_loop.h>
+#include <igl/cotmatrix.h>
+#include <igl/min_quad_with_fixed.h>
+
+using namespace std;
 
 void tutte(
   const Eigen::MatrixXd & V,
@@ -7,5 +14,55 @@ void tutte(
 {
   // Replace with your code
   U = V.leftCols(2);
-}
+  // Detect the boundary vertices
+  Eigen::VectorXi bnd;
+  igl::boundary_loop(F,bnd);
+  std::vector<std::vector<int>> lb;
+  igl::boundary_loop(F, lb);
+
+  // Map the boundary vertices to a circle
+  Eigen::MatrixXd bnd_uv;
+  igl::map_vertices_to_circle(V,bnd,bnd_uv);
+  // Changing the orientation of the boudary circle
+  bnd_uv.col(0) = -bnd_uv.col(0);
+
+  // Calculate L
+  Eigen::SparseMatrix<double> L(V.rows(), V.rows());
+  L.setZero();
+  std::vector<Eigen::Triplet<double>> triplets;
+  // Adding weights for each half edge
+  for (int f = 0; f < F.rows(); f++) {
+      int i = F(f, 1);   int j = F(f, 2);
+      double wij = 1.0/ (V.row(i) - V.row(j)).norm();
+      //wij = 1;
+      triplets.push_back({i, j,  wij });
+      triplets.push_back({j, i,  wij });
+      triplets.push_back({i, i, -wij });
+      triplets.push_back({j, j, -wij });
 
+      i = F(f, 2);   j = F(f, 0);
+      triplets.push_back({i, j,  wij });
+      triplets.push_back({j, i,  wij });
+      triplets.push_back({i, i, -wij });
+      triplets.push_back({j, j, -wij });
+
+      i = F(f, 0);   j = F(f, 1);
+      triplets.push_back({i, j,  wij });
+      triplets.push_back({j, i,  wij });
+      triplets.push_back({i, i, -wij });
+      triplets.push_back({j, j, -wij });
+  }
+  L.setFromTriplets(triplets.begin(), triplets.end());
+
+  // Calculate minimizer subject to the fixed boundary points
+  igl::min_quad_with_fixed_data<double> data;
+  Eigen::SparseMatrix<double> Aeq;
+  igl::min_quad_with_fixed_precompute(L, bnd, Aeq, false, data);
+
+  U = Eigen::MatrixXd::Zero(data.n,2);
+  Eigen::MatrixXd B(data.n, data.n);
+  B.setZero();
+  Eigen::MatrixXd Beq;
+  igl::min_quad_with_fixed_solve(data, B, bnd_uv, Beq, U);
+  igl::min_quad_with_fixed(L, B, bnd, bnd_uv, Aeq, Beq, false, U);
+}
diff --git a/src/vector_area_matrix.cpp b/src/vector_area_matrix.cpp
index 1d60ef8..a226602 100644
--- a/src/vector_area_matrix.cpp
+++ b/src/vector_area_matrix.cpp
@@ -1,4 +1,8 @@
 #include "vector_area_matrix.h"
+#include "igl/boundary_loop.h"
+#include <iostream>
+
+using namespace std;
 
 void vector_area_matrix(
   const Eigen::MatrixXi & F,
@@ -7,5 +11,27 @@ void vector_area_matrix(
   // Replace with your code
   int V_size = F.maxCoeff()+1;
   A.resize(V_size*2,V_size*2);
-}
+  A.setZero();
+  Eigen::SparseMatrix<double> A_bar(V_size * 2, V_size * 2);
+  std::vector<Eigen::Triplet<double>> triplets;
+  // Detect the boundary vertices
+  std::vector<std::vector<int>> lb;
+  igl::boundary_loop(F, lb);
 
+  for (std::vector<std::vector<int>>::iterator it = lb.begin(); it != lb.end(); ++it) {
+	  for (vector<int>::iterator bdry = (*it).begin(); bdry != (*it).end(); ++bdry) {
+		  int i, j;
+		  if ((bdry + 1) != (*it).end()) {
+			  i = *bdry; j = *(bdry+1);
+		  } else {
+			i = *bdry; j = *((*it).begin());
+		  }
+		  // Fill in entries of A_bar as well as entries of A_bar ^ T
+		  triplets.push_back({ i, j + V_size,  1./2. });
+		  triplets.push_back({ i + V_size, j, -1./2. });
+		  triplets.push_back({ j + V_size, i,  1./2. });
+		  triplets.push_back({ j, i + V_size, -1./2. });
+	  }
+  }
+  A.setFromTriplets(triplets.begin(), triplets.end());
+}