diff --git a/opty/direct_collocation.py b/opty/direct_collocation.py
index c598695e..2a236108 100644
--- a/opty/direct_collocation.py
+++ b/opty/direct_collocation.py
@@ -1181,9 +1181,6 @@ def constraints(state_values, specified_values, constant_values,
 
             """
 
-            if state_values.shape[0] < 2:
-                raise ValueError('There should always be at least two states.')
-
             assert state_values.shape == (self.num_states,
                                           self.num_collocation_nodes)
             # n x N - 1
@@ -1562,8 +1559,6 @@ def constraints_jacobian(state_values, specified_values,
             - n*(N - 1) : number of constraints
 
             """
-            if state_values.shape[0] < 2:
-                raise ValueError('There should always be at least two states.')
 
             # Each of these arrays are shape(n, N - 1). The x_adjacent is
             # either the previous value of the state or the next value of
diff --git a/opty/tests/test_direct_collocation.py b/opty/tests/test_direct_collocation.py
index adc2f3f8..f6c34b43 100644
--- a/opty/tests/test_direct_collocation.py
+++ b/opty/tests/test_direct_collocation.py
@@ -1516,7 +1516,7 @@ def test_known_and_unknown_order():
 
 def test_for_algebraic_eoms():
     """
-    If algebraic equations of motion are given to Problem, a ValueError should
+    If only algebraic equations of motion are given to Problem, a ValueError should
     be raised. This a a test for this
     """
 
@@ -1566,3 +1566,61 @@ def test_for_algebraic_eoms():
         )
 
     assert excinfo.type is ValueError
+
+def test_one_eom_only():
+    """
+    Only one differential equation should work. This tests for the corrrect
+    shape of the constraints and jacobian.
+
+    """
+    # Equations of motion.
+    t = mech.dynamicsymbols._t
+    y, u = mech.dynamicsymbols('y u')
+
+    eom = sym.Matrix([-y.diff(t) - y**3 + u])
+
+    t0, tf = 0.0, 10.0
+    num_nodes = 100
+    interval_value = (tf - t0)/(num_nodes - 1)
+
+    state_symbols = (y, )
+    specified_symbols = (u,)
+
+    # Specify the objective function and form the gradient.
+    obj_func = sym.Integral(y**2 + u**2, t)
+    obj, obj_grad = create_objective_function(
+        obj_func,
+        state_symbols,
+        specified_symbols,
+        tuple(),
+        num_nodes,
+        node_time_interval=interval_value
+    )
+
+    # Specify the symbolic instance constraints.
+    instance_constraints = (
+        y.func(t0) - 1,
+        y.func(tf) - 1.5,
+    )
+
+        # Create the optimization problem and set any options.
+    prob = Problem(
+        obj,
+        obj_grad,
+        eom,
+        state_symbols,
+        num_nodes,
+        interval_value,
+        instance_constraints=instance_constraints,
+        )
+
+    initial_guess = np.zeros(prob.num_free)
+    initial_guess[0] = 1.0
+    initial_guess[num_nodes-1] = 1.5
+
+    # assert that prob.constraints and prob.jacobian have the correct shape.
+    length = 1*(num_nodes-1) + 2
+    assert prob.constraints(initial_guess).shape == (length,)
+
+    length = (2*1 + 1 + 0 + 0) * (1*(num_nodes-1)) + 2
+    assert prob.jacobian(initial_guess).shape == (length,)