diff --git a/polytope/solvers.py b/polytope/solvers.py
index 33894ef..7cd7173 100644
--- a/polytope/solvers.py
+++ b/polytope/solvers.py
@@ -50,9 +50,17 @@
 except ImportError:
     logger.info('MOSEK solver not found.')
 
+try:
+    import gurobipy as gp
+    installed_solvers.add('gurobi')
+except ImportError:
+    logger.info('GUROBI solver not found')
+
 
 # choose default from installed choices
-if 'glpk' in installed_solvers:
+if 'gurobi' in installed_solvers:
+    default_solver = 'gurobi'
+elif 'glpk' in installed_solvers:
     default_solver = 'glpk'
 elif 'scipy' in installed_solvers:
     default_solver = 'scipy'
@@ -62,7 +70,6 @@
         "`installed_solvers` wasn't empty above?")
 
 
-
 def lpsolve(c, G, h, solver=None):
     """Try to solve linear program with given or default solver.
 
@@ -88,6 +95,8 @@ def lpsolve(c, G, h, solver=None):
         result = _solve_lp_using_cvxopt(c, G, h, solver=solver)
     elif solver == 'scipy':
         result = _solve_lp_using_scipy(c, G, h)
+    elif solver == 'gurobi':
+        result = _solve_lp_using_gurobi(c, G, h)
     else:
         raise Exception(
             'unknown LP solver "{s}".'.format(s=solver))
@@ -145,6 +154,34 @@ def _solve_lp_using_scipy(c, G, h):
         fun=sol.fun)
 
 
+def _solve_lp_using_gurobi(c, G, h):
+    """Attempt linear optimization using gurobipy"""
+    _assert_have_solver('gurobi')
+    m = gp.Model()
+    x = m.addMVar(G.shape[1], lb=-gp.GRB.INFINITY)
+    m.addConstr(G@x <= h)
+    m.setObjective(c@x)
+    m.optimize()
+
+    result = dict()
+    if m.Status == gp.GRB.OPTIMAL:
+        result['status'] = 0
+        result['x'] = x.x
+        result['fun'] = m.ObjVal
+        return result
+    elif m.Status == gp.GRB.INFEASIBLE:
+        result['status'] = 2
+    elif m.Status == gp.GRB.INF_OR_UNBD or m.Status == gp.GRB.UNBOUNDED:
+        result['status'] = 3
+    else:
+        raise ValueError(f'`gurobipy` returned unexpected status value: {m.Status}')
+
+    result['x'] = None
+    result['fun'] = None
+    return result
+
+
+
 def _assert_have_solver(solver):
     """Raise `RuntimeError` if `solver` is absent."""
     if solver in installed_solvers:
diff --git a/requirements/extras.txt b/requirements/extras.txt
index 410a6d6..d5d1d9f 100644
--- a/requirements/extras.txt
+++ b/requirements/extras.txt
@@ -1 +1 @@
-cvxopt==1.2.4
+cvxopt==1.2.4
\ No newline at end of file
diff --git a/tests/polytope_test.py b/tests/polytope_test.py
index 55c69b0..0cb0a10 100644
--- a/tests/polytope_test.py
+++ b/tests/polytope_test.py
@@ -486,5 +486,14 @@ def is_glpk_present():
         return False
 
 
+def test_gurobipy_return_same_result_as_scipy():
+    c, A, b = example_1d()
+    result_gurobi = solvers.lpsolve(c, A, b, solver='gurobi')
+    result_scipy = solvers.lpsolve(c, A, b, solver='scipy')
+    nt.assert_equal(result_scipy['status'], result_gurobi['status'])
+    nt.assert_almost_equal(result_scipy['x'][0], result_gurobi['x'][0])
+    nt.assert_almost_equal(result_scipy['fun'], result_gurobi['fun'])
+
+
 if __name__ == '__main__':
     pass