diff --git a/include/igl/ismember.cpp b/include/igl/ismember.cpp
new file mode 100644
index 000000000..3a3f927eb
--- /dev/null
+++ b/include/igl/ismember.cpp
@@ -0,0 +1,96 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "ismember.h"
+#include "colon.h"
+#include "list_to_matrix.h"
+#include "sort.h"
+#include "unique.h"
+
+template <
+  typename DerivedA,
+  typename DerivedB,
+  typename DerivedIA,
+  typename DerivedLOCB>
+IGL_INLINE void igl::ismember(
+  const Eigen::PlainObjectBase<DerivedA> & A,
+  const Eigen::PlainObjectBase<DerivedB> & B,
+  Eigen::PlainObjectBase<DerivedIA> & IA,
+  Eigen::PlainObjectBase<DerivedLOCB> & LOCB)
+{
+  using namespace Eigen;
+  using namespace std;
+  IA.resize(A.rows(),A.cols());
+  IA.setConstant(false);
+  LOCB.resize(A.rows(),A.cols());
+  LOCB.setConstant(-1);
+  // boring base cases
+  if(A.size() == 0)
+  {
+    return;
+  }
+  if(B.size() == 0)
+  {
+    return;
+  }
+
+  // Get rid of any duplicates
+  typedef Matrix<typename DerivedA::Scalar,Dynamic,1> VectorA;
+  typedef Matrix<typename DerivedB::Scalar,Dynamic,1> VectorB;
+  const VectorA vA(Eigen::Map<const VectorA>(A.data(), A.cols()*A.rows(),1));
+  const VectorB vB(Eigen::Map<const VectorB>(B.data(), B.cols()*B.rows(),1));
+  VectorA uA;
+  VectorB uB;
+  Eigen::Matrix<typename DerivedA::Index,Dynamic,1> uIA,uIuA,uIB,uIuB;
+  unique(vA,uA,uIA,uIuA);
+  unique(vB,uB,uIB,uIuB);
+  // Sort both
+  VectorA sA;
+  VectorB sB;
+  Eigen::Matrix<typename DerivedA::Index,Dynamic,1> sIA,sIB;
+  sort(uA,1,true,sA,sIA);
+  sort(uB,1,true,sB,sIB);
+
+  Eigen::Matrix<bool,Eigen::Dynamic,1> uF = 
+    Eigen::Matrix<bool,Eigen::Dynamic,1>::Zero(sA.size(),1);
+  Eigen::Matrix<typename DerivedLOCB::Scalar, Eigen::Dynamic,1> uLOCB =
+    Eigen::Matrix<typename DerivedLOCB::Scalar,Eigen::Dynamic,1>::
+    Constant(sA.size(),1,-1);
+  {
+    int bi = 0;
+    // loop over sA
+    bool past = false;
+    for(int a = 0;a<sA.size();a++)
+    {
+      while(!past && sA(a)>sB(bi))
+      {
+        bi++;
+        past = bi>=sB.size();
+      }
+      if(!past && sA(a)==sB(bi))
+      {
+        uF(sIA(a)) = true;
+        uLOCB(sIA(a)) = uIB(sIB(bi));
+      }
+    }
+  }
+
+  Map< Matrix<typename DerivedIA::Scalar,Dynamic,1> > 
+    vIA(IA.data(),IA.cols()*IA.rows(),1);
+  Map< Matrix<typename DerivedLOCB::Scalar,Dynamic,1> > 
+    vLOCB(LOCB.data(),LOCB.cols()*LOCB.rows(),1);
+  for(int a = 0;a<A.size();a++)
+  {
+    vIA(a) = uF(uIuA(a));
+    vLOCB(a) = uLOCB(uIuA(a));
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::ismember<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif
diff --git a/include/igl/ismember.h b/include/igl/ismember.h
new file mode 100644
index 000000000..7c1889dae
--- /dev/null
+++ b/include/igl/ismember.h
@@ -0,0 +1,41 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_ISMEMBER_H
+#define IGL_ISMEMBER_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Determine if elements of A exist in elements of B
+  //
+  // Inputs:
+  //   A  ma by na matrix
+  //   B  mb by nb matrix
+  // Outputs:
+  //   IA  ma by na matrix of flags whether corresponding element of A exists in
+  //     B
+  //   LOCB  ma by na matrix of indices in B locating matching element (-1 if
+  //     not found), indices assume column major ordering
+  //
+  template <
+    typename DerivedA,
+    typename DerivedB,
+    typename DerivedIA,
+    typename DerivedLOCB>
+  IGL_INLINE void ismember(
+    const Eigen::PlainObjectBase<DerivedA> & A,
+    const Eigen::PlainObjectBase<DerivedB> & B,
+    Eigen::PlainObjectBase<DerivedIA> & IA,
+    Eigen::PlainObjectBase<DerivedLOCB> & LOCB);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "ismember.cpp"
+#endif
+#endif
+
diff --git a/matlab-to-eigen.html b/matlab-to-eigen.html
index c8b83f3cb..54e135b1b 100644
--- a/matlab-to-eigen.html
+++ b/matlab-to-eigen.html
@@ -272,6 +272,12 @@ bool a = A.array().any();
         <td></td>
       </tr>
 
+      <tr class=d1>
+        <td><pre><code>[IA,LOCB] = ismember(A,B)</code></pre></td>
+        <td><pre><code>igl::ismember(A,B,IA,LOCB)</code></pre></td>
+        <td></td>
+      </tr>
+
       <!-- Insert rows for each command pair -->
 
       <!-- Leave this here for copy and pasting
diff --git a/tutorial/tutorial.md b/tutorial/tutorial.md
index 86fc9b45d..7e9a08456 100644
--- a/tutorial/tutorial.md
+++ b/tutorial/tutorial.md
@@ -877,6 +877,7 @@ functionality as common Matlab functions.
 | `igl::hsv_to_rgb`        | Convert HSV colors to RGB (cf. Matlab's `hsv2rgb`) |
 | `igl::intersect`         | Set intersection of matrix elements. |
 | `igl::isdiag`            | Determine whether matrix is diagonal |
+| `igl::ismember`          | Determine whether elements in A occur in B |
 | `igl::jet`               | Quantized colors along the rainbow. |
 | `igl::median`            | Compute the median per column |
 | `igl::mode`              | Compute the mode per column |