diff --git a/pybamm/solvers/c_solvers/idaklu/casadi_solver.cpp b/pybamm/solvers/c_solvers/idaklu/casadi_solver.cpp
index e30a59d852..3cb72b67a6 100644
--- a/pybamm/solvers/c_solvers/idaklu/casadi_solver.cpp
+++ b/pybamm/solvers/c_solvers/idaklu/casadi_solver.cpp
@@ -344,25 +344,12 @@ Solution CasadiSolver::solve(np_array t_np, np_array y0_np, np_array yp0_np,
     IDASensReInit(ida_mem, IDA_SIMULTANEOUS, yyS, ypS);
   }
 
-  // TODO: not sure if I need to do this again...
-  for (int is = 0 ; is < number_of_parameters; is++) {
-    ySval[is] = N_VGetArrayPointer(yyS[is]);
-  }
-  yval = N_VGetArrayPointer(yy);
-  ypval = N_VGetArrayPointer(yp);
-
-
   // calculate consistent initial conditions
   DEBUG("IDACalcIC");
   IDACalcIC(ida_mem, IDA_YA_YDP_INIT, t(1));
-  for (int j = 0; j < number_of_parameters; j++) {
-    for (int k = 0; k < number_of_states; k++) {
-      std::cout << "ySval[" << j << "]["<< k << "] = " << ySval[j][k] << std::endl;
-    }
-    //DEBUG_VECTOR(yyS[j]);
-  }
-  for (int k = 0; k < 10; k++) {
-    std::cout << "yval["<< k << "] = " << yval[k] << std::endl;
+  if (number_of_parameters > 0)
+  {
+    IDAGetSens(ida_mem, &t0, yyS);
   }
 
   int t_i = 1;
diff --git a/pybamm/solvers/c_solvers/idaklu/casadi_sundials_functions.cpp b/pybamm/solvers/c_solvers/idaklu/casadi_sundials_functions.cpp
index cd54ad15e7..78a5138342 100644
--- a/pybamm/solvers/c_solvers/idaklu/casadi_sundials_functions.cpp
+++ b/pybamm/solvers/c_solvers/idaklu/casadi_sundials_functions.cpp
@@ -279,7 +279,7 @@ int sensitivities_casadi(int Ns, realtype t, N_Vector yy, N_Vector yp,
   }
   // resvalsS now has (∂F/∂p i )
   p_python_functions->sens();
-
+  
   for (int i = 0; i < np; i++)
   {
     // put (∂F/∂y)s i (t) in tmp
diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py
index b26ce8262d..743b04f8ec 100644
--- a/pybamm/solvers/idaklu_solver.py
+++ b/pybamm/solvers/idaklu_solver.py
@@ -501,7 +501,6 @@ def _integrate(self, model, t_eval, inputs_dict=None):
 
             y0S = (x.full() for x in y0S)
             y0S = [x.flatten() for x in y0S]
-        print("integrate y0S", y0S)
 
         # solver works with ydot0 set to zero
         ydot0 = np.zeros_like(y0)
@@ -523,7 +522,6 @@ def _integrate(self, model, t_eval, inputs_dict=None):
 
         timer = pybamm.Timer()
         if model.convert_to_format == "casadi":
-            print("sol for ", t_eval)
             sol = self._setup["solver"].solve(
                 t_eval,
                 y0full,
diff --git a/tests/integration/test_models/standard_model_tests.py b/tests/integration/test_models/standard_model_tests.py
index 3c1b6eba24..9341122d84 100644
--- a/tests/integration/test_models/standard_model_tests.py
+++ b/tests/integration/test_models/standard_model_tests.py
@@ -122,7 +122,7 @@ def test_sensitivities(self, param_name, param_value, output_name="Voltage [V]")
         output_sens = self.solution[output_name].sensitivities[param_name]
 
         # check via finite differencing
-        h = 1e-4 * param_value
+        h = 1e-2 * param_value
         inputs_plus = {param_name: (param_value + 0.5 * h)}
         inputs_neg = {param_name: (param_value - 0.5 * h)}
         sol_plus = self.solver.solve(self.model, t_eval, inputs=inputs_plus)
diff --git a/tests/unit/test_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py
index 4048a437f7..efd0439f32 100644
--- a/tests/unit/test_solvers/test_idaklu_solver.py
+++ b/tests/unit/test_solvers/test_idaklu_solver.py
@@ -204,10 +204,8 @@ def test_sensitivites_initial_condition(self):
             model, t_eval, inputs={"a": a_value}, calculate_sensitivities=True
         )
 
-        print(sol["2v"].sensitivities["a"].full().flatten())
-
         np.testing.assert_array_almost_equal(
-            sol["2v"].sensitivities["a"].full().flatten(), np.exp(-sol.t) * 2
+            sol["2v"].sensitivities["a"].full().flatten(), np.exp(-sol.t) * 2, decimal=4
         )
 
     def test_ida_roberts_klu_sensitivities(self):