diff --git a/CHANGELOG.md b/CHANGELOG.md
index 42314363b1..4fbc229020 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -20,7 +20,7 @@
 - Removed the `start_step_offset` setting and disabled minimum `dt` warnings for drive cycles with the (`IDAKLUSolver`). ([#4416](https://github.com/pybamm-team/PyBaMM/pull/4416))
 
 ## Bug Fixes
-
+- Fixed bug in post-processing solutions with infeasible experiments using the (`IDAKLUSolver`). ([#4541](https://github.com/pybamm-team/PyBaMM/pull/4541))
 - Disabled IREE on MacOS due to compatibility issues and added the CasADI
   path to the environment to resolve issues on MacOS and Linux. Windows
   users may still experience issues with interpolation. ([#4528](https://github.com/pybamm-team/PyBaMM/pull/4528))
diff --git a/src/pybamm/solvers/c_solvers/idaklu/observe.cpp b/src/pybamm/solvers/c_solvers/idaklu/observe.cpp
index 8f1d90e55d..6661bdecc9 100644
--- a/src/pybamm/solvers/c_solvers/idaklu/observe.cpp
+++ b/src/pybamm/solvers/c_solvers/idaklu/observe.cpp
@@ -100,6 +100,11 @@ class TimeSeriesInterpolator {
         ) {
         for (size_t i = 0; i < ts_data_np.size(); i++) {
             const auto& t_data = ts_data_np[i].unchecked<1>();
+            // Continue if there is no data
+            if (t_data.size() == 0) {
+                continue;
+            }
+
             const realtype t_data_final = t_data(t_data.size() - 1);
             realtype t_interp_next = t_interp(i_interp);
             // Continue if the next interpolation point is beyond the final data point
@@ -227,6 +232,10 @@ class TimeSeriesProcessor {
         int i_entries = 0;
         for (size_t i = 0; i < ts.size(); i++) {
             const auto& t = ts[i].unchecked<1>();
+            // Continue if there is no data
+            if (t.size() == 0) {
+                continue;
+            }
             const auto& y = ys[i].unchecked<2>();
             const auto input = inputs[i].data();
             const auto func = *funcs[i];
diff --git a/src/pybamm/solvers/processed_variable.py b/src/pybamm/solvers/processed_variable.py
index 12cf140b38..3de6e4bd50 100644
--- a/src/pybamm/solvers/processed_variable.py
+++ b/src/pybamm/solvers/processed_variable.py
@@ -133,19 +133,22 @@ def _setup_cpp_inputs(self, t, full_range):
         ys = self.all_ys
         yps = self.all_yps
         inputs = self.all_inputs_casadi
-        # Find the indices of the time points to observe
-        if full_range:
-            idxs = range(len(ts))
-        else:
-            idxs = _find_ts_indices(ts, t)
 
-        if isinstance(idxs, list):
-            # Extract the time points and inputs
-            ts = [ts[idx] for idx in idxs]
-            ys = [ys[idx] for idx in idxs]
-            if self.hermite_interpolation:
-                yps = [yps[idx] for idx in idxs]
-            inputs = [self.all_inputs_casadi[idx] for idx in idxs]
+        # Remove all empty ts
+        idxs = np.where([ti.size > 0 for ti in ts])[0]
+
+        # Find the indices of the time points to observe
+        if not full_range:
+            ts_nonempty = [ts[idx] for idx in idxs]
+            idxs_subset = _find_ts_indices(ts_nonempty, t)
+            idxs = idxs[idxs_subset]
+
+        # Extract the time points and inputs
+        ts = [ts[idx] for idx in idxs]
+        ys = [ys[idx] for idx in idxs]
+        if self.hermite_interpolation:
+            yps = [yps[idx] for idx in idxs]
+        inputs = [self.all_inputs_casadi[idx] for idx in idxs]
 
         is_f_contiguous = _is_f_contiguous(ys)
 
@@ -977,8 +980,4 @@ def _find_ts_indices(ts, t):
     if (t[-1] > ts[-1][-1]) and (len(indices) == 0 or indices[-1] != len(ts) - 1):
         indices.append(len(ts) - 1)
 
-    if len(indices) == len(ts):
-        # All indices are included
-        return range(len(ts))
-
     return indices
diff --git a/src/pybamm/solvers/solution.py b/src/pybamm/solvers/solution.py
index 1aa540ab3c..7410c2161a 100644
--- a/src/pybamm/solvers/solution.py
+++ b/src/pybamm/solvers/solution.py
@@ -580,11 +580,15 @@ def _update_variable(self, variable):
         # Iterate through all models, some may be in the list several times and
         # therefore only get set up once
         vars_casadi = []
-        for i, (model, ys, inputs, var_pybamm) in enumerate(
-            zip(self.all_models, self.all_ys, self.all_inputs, vars_pybamm)
+        for i, (model, ts, ys, inputs, var_pybamm) in enumerate(
+            zip(self.all_models, self.all_ts, self.all_ys, self.all_inputs, vars_pybamm)
         ):
-            if ys.size == 0 and var_pybamm.has_symbol_of_classes(
-                pybamm.expression_tree.state_vector.StateVector
+            if (
+                ys.size == 0
+                and var_pybamm.has_symbol_of_classes(
+                    pybamm.expression_tree.state_vector.StateVector
+                )
+                and not ts.size == 0
             ):
                 raise KeyError(
                     f"Cannot process variable '{variable}' as it was not part of the "