diff --git a/symfit/core/models.py b/symfit/core/models.py
index 5563ebf6..4d373cc3 100644
--- a/symfit/core/models.py
+++ b/symfit/core/models.py
@@ -1167,17 +1167,17 @@ def eval_components(self, *args, **kwargs):
         # np.flip is needed because odeint wants the first point to be t_initial
         # and so t_smaller is a declining series.
         if t_initial in t_like:
-            t_bigger = t_like[t_like >= t_initial]
+            t_bigger = t_like[t_like > t_initial]
             t_smaller = t_like[t_like <= t_initial][::-1]
         else:
             t_bigger = np.concatenate(
                 (np.array([t_initial]), t_like[t_like > t_initial])
             )
             t_smaller = np.concatenate(
-                (np.array([t_initial]), t_like[t_like < t_initial][::-1])
+                (np.array([t_initial]), t_like[t_like <= t_initial][::-1])
             )
         # Properly ordered time axis containing t_initial
-        t_total = np.concatenate((t_smaller[::-1][:-1], t_bigger))
+        t_total = np.concatenate((t_smaller[::-1], t_bigger))
 
         # Call the numerical integrator. Note that we only pass the
         # model_params, which will be used by sympy_to_py to create something we
@@ -1203,7 +1203,7 @@ def eval_components(self, *args, **kwargs):
             *self.lsoda_args, **self.lsoda_kwargs
         )
 
-        ans = np.concatenate((ans_smaller[1:][::-1], ans_bigger))
+        ans = np.concatenate((ans_smaller[::-1], ans_bigger))
         if t_initial in t_like:
             # The user also requested to know the value at t_initial, so keep it.
             return ans.T
diff --git a/tests/test_ode.py b/tests/test_ode.py
index f53b2041..6e510065 100644
--- a/tests/test_ode.py
+++ b/tests/test_ode.py
@@ -258,3 +258,27 @@ def test_initial_parameters():
     assert ode_model.params == [a0, c0, k, l, m, p]
     assert ode_model.initial_params == [a0, c0]
     assert ode_model.model_params == [a0, k, l, m, p]
+
+
+def test_duplicate_initial_t():
+    """
+    Tests that having duplicate t0 values in tdata doesn't cause shape mismatch issues, such as in #339
+    """
+    tdata = np.array([0, 10, 10, 26, 44, 70, 120])
+    adata = 10e-4 * np.array([54, 54, 44, 34, 27, 20, 14])
+
+    a, b, t = variables('a, b, t')
+    k = Parameter('k', 0.1)
+
+    model_dict = {
+        D(a, t): - k * a ** 2,
+        D(b, t): k * a ** 2,
+    }
+
+    ode_model = ODEModel(model_dict, initial={t: 26, a: adata[0], b: 0.0})
+
+    fit = Fit(ode_model, t=tdata, a=adata, b=None)
+    try:
+        fit_result = fit.execute()
+    except ValueError:
+        pytest.fail()