diff --git a/pyoptmat/chunktime.py b/pyoptmat/chunktime.py
index 6187bc0..f9f3ec4 100644
--- a/pyoptmat/chunktime.py
+++ b/pyoptmat/chunktime.py
@@ -50,6 +50,7 @@ def newton_raphson_chunk(
     i = 0
 
     while (i < miter) and torch.any(torch.logical_not(torch.logical_or(nR <= atol, nR / nR0 <= rtol))):
+        print(i, torch.max(nR))
         dx = solver.solve(J, R)
         if linesearch:
             x, R, J, nR, alpha = chunk_linesearch(x, dx, fn, R, rtol, atol)
diff --git a/pyoptmat/flowrules.py b/pyoptmat/flowrules.py
index c9034aa..595d615 100644
--- a/pyoptmat/flowrules.py
+++ b/pyoptmat/flowrules.py
@@ -40,9 +40,6 @@ class FlowRule(nn.Module):
     def __init__(self):
         super().__init__()
 
-    def setup(self, t, y):
-        pass
-
     def dflow_dhist(self, s, h, t, T, e):
         """
         The derivative of the flow rate with respect to the internal variables
@@ -190,10 +187,6 @@ def __init__(
                 "model must have the same number of internal variables"
             )
 
-    def setup(self, t, y):
-        self.model1.setup(t, y)
-        self.model2.setup(t, y)
-
     @property
     def nhist(self):
         """
@@ -459,10 +452,6 @@ def __init__(
                 "model must have the same number of internal variables"
             )
 
-    def setup(self, t, y):
-        self.model1.setup(t, y)
-        self.model2.setup(t, y)
-
     @property
     def nhist(self):
         """
@@ -1541,15 +1530,25 @@ def dhist_derate(self, s, h, t, T, e):
 
 class IsoKinRateIndependentPlasticity(FlowRule):
     """
-    Viscoplasticity with isotropic and kinematic hardening, defined as
+    An approximation to rate independent plasticity.  The model is defined by 
+    
+    .. math::
+
+        \\dot{\\varepsilon}_{in} = \\xi(f) \\dot{\\varepsilon}_{p,ri}
+
+    where :math:`\\xi` is a sigmoid function, :math:`f` is a flow surface
+    and :math:`\\dot{\\varepsilon}_{p,ri}` is the rate independent plastic flow
+    rate, as defined by the classical consistency conditions.
+
+    This function uses the flow surface
 
     .. math::
 
-      \\dot{\\varepsilon}_{in}=\\left\\langle \\frac{\\left|\\sigma-x\\right|-s_{0}-k}{\\eta}\\right\\rangle ^{n}\\operatorname{sign}\\left(\\sigma-X\\right)
+        f = \\left|\\sigma - x\\right| - \\sigma_y - k
 
     and where the :py:class:`pyoptmat.hardening.IsotropicHardeningModel` and
     :py:class:`pyoptmat.hardening.KinematicHardeningModel` objects determine the
-    history rate.
+    history rate of isotropic (:math:`k`) and kinematic (:math:`x`) hardening.
 
     The :py:class:`pyoptmat.hardening.IsotropicHardeningModel` and
     :py:class:`pyoptmat.hardening.KinematicHardeningModel` objects each define both a
@@ -1562,14 +1561,17 @@ class IsoKinRateIndependentPlasticity(FlowRule):
     the information this class needs to provide.
 
     Args:
-      n (|TP|):         rate sensitivity
-      eta (|TP|):       flow viscosity
-      s0 (|TP|):        initial value of flow stress (i.e. the threshold stress)
+      E (|TP|):         young's modulus, needed to define the rate independent flow rate
+      sy (|TP|):        yield stress
       isotropic (:py:class:`pyoptmat.hardening.IsotropicHardeningModel`): object providing the isotropic hardening model
       kinematic (:py:class:`pyoptmat.hardening.IsotropicHardeningModel`): object providing the kinematic hardening model
+
+    Keyword Args:
+        soffset (float):    small offset to the stress in the  yield surface to avoid a singularity at zero
+        s (float):          scale factor for the sigmoid function, controls the amount of smoothing at the onset of plasticity
     """
 
-    def __init__(self, E, sy, isotropic, kinematic, soffset = 1e-10):
+    def __init__(self, E, sy, isotropic, kinematic, soffset = 1e-10, s = 1.0):
         super().__init__()
         self. E = E
         self.isotropic = isotropic
@@ -1577,16 +1579,7 @@ def __init__(self, E, sy, isotropic, kinematic, soffset = 1e-10):
         self.sy = sy
         self.soffset = soffset
 
-        self.yielding = None
-
-    def setup(self, t, y):
-        self.yielding = None 
-
-    def update_yielding(self, s, h, T):
-        if self.yielding is None or self.yielding.shape != s.shape:
-            self.yielding = self.f(s, h, T) > 0
-        else:
-            self.yielding = torch.logical_or(self.yielding,self.f(s, h, T) > 0)
+        self.s = s
 
     def f(self, s, h, T):
         """
@@ -1598,16 +1591,60 @@ def f(self, s, h, T):
           T (torch.tensor):   temperature
 
         Returns:
-          the value of the yield function...
+          the value of the yield function
         """
         ih = self.isotropic.value(h[..., : self.isotropic.nhist])
         kh = self.kinematic.value(h[..., self.isotropic.nhist :])
         return torch.abs(s-kh) - self.sy(T) - ih
 
+    def df_ds(self, s, h, T):
+        """
+        The derivative of the yield function with respect to stress
+
+        Args:
+          s (torch.tensor):   stress
+          h (torch.tensor):   history
+          T (torch.tensor):   temperature
+
+        Returns:
+          the derivative value
+        """
+        kh = self.kinematic.value(h[..., self.isotropic.nhist :])
+        return torch.sign(s-kh)
+
+    def df_dh(self, s, h, T):
+        """
+        The derivative of the yield function with respect to the history
+
+        Args:
+          s (torch.tensor):   stress
+          h (torch.tensor):   history
+          T (torch.tensor):   temperature
+
+        Returns:
+          the derivative value
+        """
+        kh = self.kinematic.value(h[..., self.isotropic.nhist :])
+        di = -self.isotropic.dvalue(h[..., : self.isotropic.nhist])
+        dk = -torch.sign(s-kh).unsqueeze(-1) * self.kinematic.dvalue(h[..., self.isotropic.nhist :])
+
+        return torch.cat([di,dk], axis = -1)
+
     def ep_residual(self, ep, s, h, t, T, e):
         """
         The residual function to solve for the consistency parameter, here just
         expanded to the plastic flow rate as it's a scalar.
+
+        Args:
+          ep (torch.tensor):  current values of the plastic flow rate
+          s (torch.tensor):   stress
+          h (torch.tensor):   history
+          t (torch.tensor):   current time
+          T (torch.tensor):   temperature
+          e (torch.tensor):   current total strain rate
+
+        Returns:
+          the residual value
         """
         hiso = h[..., : self.isotropic.nhist]
         hkin = h[..., self.isotropic.nhist :]
@@ -1621,7 +1658,18 @@ def ep_residual(self, ep, s, h, t, T, e):
 
     def ep_jacobian_ep(self, ep, s, h, t, T, e):
         """
-        The Jacobian for the plastic flow rate residual with respect to the plastic strain rate
+        The Jacobian of the consistency residual with respect to the plastic strain rate.
+
+        Args:
+          ep (torch.tensor):  current values of the plastic flow rate
+          s (torch.tensor):   stress
+          h (torch.tensor):   history
+          t (torch.tensor):   current time
+          T (torch.tensor):   temperature
+          e (torch.tensor):   current total strain rate
+
+        Returns:
+          the derivative value
         """
         hiso = h[..., : self.isotropic.nhist]
         hkin = h[..., self.isotropic.nhist :]
@@ -1635,7 +1683,18 @@ def ep_jacobian_ep(self, ep, s, h, t, T, e):
 
     def ep_jacobian_s(self, ep, s, h, t, T, e):
         """
-        The Jacobian for the plastic flow rate residual with respect to the stress
+        The Jacobian of the consistency residual with respect to the stress
+
+        Args:
+          ep (torch.tensor):  current values of the plastic flow rate
+          s (torch.tensor):   stress
+          h (torch.tensor):   history
+          t (torch.tensor):   current time
+          T (torch.tensor):   temperature
+          e (torch.tensor):   current total strain rate
+
+        Returns:
+          the derivative value
         """
         hiso = h[..., : self.isotropic.nhist]
         hkin = h[..., self.isotropic.nhist :]
@@ -1649,7 +1708,18 @@ def ep_jacobian_s(self, ep, s, h, t, T, e):
 
     def ep_jacobian_h(self, ep, s, h, t, T, e):
         """
-        The Jacobian for the plastic flow rate residual with respect to the stress
+        The Jacobian of the consistency residual with respect to the history
+
+        Args:
+          ep (torch.tensor):  current values of the plastic flow rate
+          s (torch.tensor):   stress
+          h (torch.tensor):   history
+          t (torch.tensor):   current time
+          T (torch.tensor):   temperature
+          e (torch.tensor):   current total strain rate
+
+        Returns:
+          the derivative value
         """
         hiso = h[..., : self.isotropic.nhist]
         hkin = h[..., self.isotropic.nhist :]
@@ -1663,7 +1733,18 @@ def ep_jacobian_h(self, ep, s, h, t, T, e):
 
     def ep_jacobian_e(self, ep, s, h, t, T, e):
         """
-        The Jacobian for the plastic flow rate residual with respect to the total strain rate
+        The Jacobian of the consistency residual with respect to the total strain rate
+
+        Args:
+          ep (torch.tensor):  current values of the plastic flow rate
+          s (torch.tensor):   stress
+          h (torch.tensor):   history
+          t (torch.tensor):   current time
+          T (torch.tensor):   temperature
+          e (torch.tensor):   current total strain rate
+
+        Returns:
+          the derivative value
         """
         hkin = h[..., self.isotropic.nhist :]
 
@@ -1671,6 +1752,79 @@ def ep_jacobian_e(self, ep, s, h, t, T, e):
         
         return (torch.sign(s - kh + self.soffset) * self.E(T)).unsqueeze(-1)
 
+    def sig(self, x):
+        """
+        Our chosen sigmoid function
+
+        .. math::
+
+            \\xi = \\frac{\\tanh(s x) + 1}{2}
+
+        with :math:`s` the scale factor
+
+        Args:
+            x (torch.tensor): input to sigmoid
+        """
+        return (torch.tanh(self.s * x) + 1.0)/2.0
+
+    def dsig(self, x):
+        """
+        Derivative of the sigmoid
+
+        Args:
+            x (torch.tensor): input to sigmoid
+        """
+        return 0.5*self.s / torch.cosh(self.s * x)**2.0
+
+    def mix_fn(self, s, h, T):
+        """
+        The mixture function :math:`\\xi(f)`
+
+        Args:
+          s (torch.tensor):   stress
+          h (torch.tensor):   history
+          T (torch.tensor):   temperature
+        """
+        return self.sig(self.f(s, h, T))
+
+    def dmix_fn_ds(self, s, h, T):
+        """
+        The derivative of the mixture function with respect to the stress
+
+        Args:
+          s (torch.tensor):   stress
+          h (torch.tensor):   history
+          T (torch.tensor):   temperature
+        """
+        return self.dsig(self.f(s, h, T)) * self.df_ds(s, h, T)
+
+    def dmix_fn_dh(self, s, h, T):
+        """
+        The derivative of the mixture function with respect to the history
+
+        Args:
+          s (torch.tensor):   stress
+          h (torch.tensor):   history
+          T (torch.tensor):   temperature
+        """
+        return self.dsig(self.f(s, h, T)).unsqueeze(-1) * self.df_dh(s, h, T)
+
+    def plastic_rate(self, s, h, t, T, e):
+        """
+        Solve for the plastic strain rate that meets the consistency criteria
+
+        Args:
+          s (torch.tensor):   stress
+          h (torch.tensor):   history
+          t (torch.tensor):   time
+          T (torch.tensor):   temperature
+          e (torch.tensor):   total strain rate
+        """
+        def RJ(x):
+            return self.ep_residual(x, s+self.soffset, h, t, T, e), self.ep_jacobian_ep(x, s+self.soffset, h, t, T, e).squeeze(-1)
+
+        return solvers.scalar_newton(RJ, torch.zeros_like(e))
+
     def flow_rate(self, s, h, t, T, e):
         """
         The flow rate itself and the derivative with respect to stress
@@ -1686,17 +1840,13 @@ def flow_rate(self, s, h, t, T, e):
           tuple(torch.tensor, torch.tensor):  the flow rate and the derivative of the flow rate with
                                               respect to stress
         """
-        fr = torch.zeros_like(e)
-        self.update_yielding(s, h, T)
-       
-        def RJ(x):
-            return self.ep_residual(x, s+self.soffset, h, t, T, e), self.ep_jacobian_ep(x, s+self.soffset, h, t, T, e).squeeze(-1)
-
-        pfr = solvers.scalar_newton(RJ, torch.zeros_like(e))
-        fr[self.yielding] = pfr[self.yielding]
-        dfr = torch.zeros_like(fr)
-        dfr[self.yielding] = -(self.ep_jacobian_s(pfr, s, h, t, T, e) / self.ep_jacobian_ep(pfr, s, h, t, T, e)).squeeze(-1)[self.yielding]
+        # Solve for the plastic flow rate
+        mf = self.mix_fn(s, h, T)
+        bfr = self.plastic_rate(s, h, t, T, e)
 
+        fr = mf * bfr
+        dfr = -(self.ep_jacobian_s(bfr, s, h, t, T, e) / self.ep_jacobian_ep(bfr, s, h, t, T, e)).squeeze(-1) * mf + self.dmix_fn_ds(s, h, T) * bfr
+        
         return fr, dfr
         
 
@@ -1718,9 +1868,9 @@ def dflow_derate(self, s, h, t, T, e):
           torch.tensor:       derivative of flow rate with respect to the
                               internal variables
         """
-        pfr = self.flow_rate(s, h, t, T, e)[0]
-        dfr = torch.zeros_like(e)
-        dfr[self.yielding] = -(self.ep_jacobian_e(pfr, s, h, t, T, e) / self.ep_jacobian_ep(pfr, s, h, t, T, e)).squeeze(-1)[self.yielding]
+        pfr = self.plastic_rate(s, h, t, T, e)
+        mf = self.mix_fn(s, h, T)
+        dfr = -(self.ep_jacobian_e(pfr, s, h, t, T, e) / self.ep_jacobian_ep(pfr, s, h, t, T, e)).squeeze(-1) * mf
 
         return dfr
 
@@ -1738,11 +1888,10 @@ def dflow_dhist(self, s, h, t, T, e):
         Returns:
           torch.tensor:       the derivative of the flow rate
         """
-        pfr = self.flow_rate(s, h, t, T, e)[0]
-        dfr = torch.zeros_like(h)
-
-        dfr[self.yielding] = -(self.ep_jacobian_h(pfr, s, h, t, T, e) / self.ep_jacobian_ep(pfr, s, h, t, T, e))[self.yielding] 
-
+        bfr = self.plastic_rate(s, h, t, T, e)
+        mf = self.mix_fn(s, h, T)
+        dfr = -(self.ep_jacobian_h(bfr, s, h, t, T, e) / self.ep_jacobian_ep(bfr, s, h, t, T, e)) * mf.unsqueeze(-1) + self.dmix_fn_dh(s, h, T) * bfr.unsqueeze(-1)
+        
         return dfr.unsqueeze(-2)
 
     @property
diff --git a/pyoptmat/models.py b/pyoptmat/models.py
index 93cdfe1..09f6660 100644
--- a/pyoptmat/models.py
+++ b/pyoptmat/models.py
@@ -72,9 +72,6 @@ def nhist(self):
         """
         return 1 + self.flowrule.nhist
 
-    def setup(self, t, y):
-        self.flowrule.setup(t, y)
-
     def forward(self, t, y, erate, T):
         """
         Return the rate equations for the strain-based version of the model
@@ -91,8 +88,8 @@ def forward(self, t, y, erate, T):
           d_y_dot_d_erate:(nbatch,1+nhist) Jacobian wrt the strain rate
           d_y_dot_d_T:    (nbatch,1+nhist) derivative wrt temperature (unused)
         """
-        stress = y[..., 0].clone()
-        h = y[..., 1:].clone()
+        stress = y[..., 0]
+        h = y[..., 1:]
 
         frate, dfrate = self.flowrule.flow_rate(stress, h, t, T, erate)
         hrate, dhrate = self.flowrule.history_rate(stress, h, t, T, erate)
@@ -455,12 +452,6 @@ def __init__(self, model, times, rates, base, temps, control, bisect_first = Fal
              bisect_first = bisect_first
         )
 
-    def setup(self, t, y):
-        if torch.any(self.econtrol):
-            self.emodel.setup(t[...,self.econtrol], y[...,self.econtrol,:])
-        if torch.any(self.scontrol):
-            self.smodel.setup(t[...,self.scontrol], y[...,self.scontrol,:])
-
     def forward(self, t, y):
         """
         Evaluate both strain and stress control and paste into the right
@@ -505,9 +496,6 @@ def __init__(self, model, erate_fn, T_fn, *args, **kwargs):
         self.erate_fn = erate_fn
         self.T_fn = T_fn
 
-    def setup(self, t, y):
-        self.model.setup(t, y)
-
     def forward(self, t, y):
         """
         Strain rate as a function of t and state
@@ -541,9 +529,6 @@ def __init__(self, model, srate_fn, stress_fn, T_fn, min_erate = -1e2, max_erate
         self.max_erate = max_erate
         self.bisect_first = bisect_first
 
-    def setup(self, t, y):
-        self.model.setup(t, y)
-
     def forward(self, t, y):
         """
         Stress rate as a function of t and state
@@ -559,7 +544,7 @@ def forward(self, t, y):
         def RJ(erate):
             ydot, _, Je, _ = self.model(t, 
                     torch.cat([cs.unsqueeze(-1), y[...,1:]], dim = -1),
-                    erate.detach(), cT)
+                    erate, cT)
 
             R = ydot[..., 0] - csr
             J = Je[..., 0]
diff --git a/pyoptmat/ode.py b/pyoptmat/ode.py
index 716cf27..ba8940b 100644
--- a/pyoptmat/ode.py
+++ b/pyoptmat/ode.py
@@ -521,10 +521,6 @@ def block_update(self, t, t_start, y_start, func, y_guess):
         # Various useful sizes
         n = t.shape[0]  # Number of time steps to do at once
         
-        self.func.setup(torch.vstack((t_start.unsqueeze(0), t)),
-                torch.vstack((torch.zeros_like(y_start).unsqueeze(0), y_guess.reshape(self.batch_size, n, self.prob_size).transpose(0, 1)))
-                + y_start.unsqueeze(0).expand(n + 1, -1, -1))
-
         def RJ(dy):
             # Make things into a more rational shape
             dy = dy.reshape(self.batch_size, n, self.prob_size).transpose(0, 1)
diff --git a/pyoptmat/solvers.py b/pyoptmat/solvers.py
index bcaf7d4..8deed4a 100644
--- a/pyoptmat/solvers.py
+++ b/pyoptmat/solvers.py
@@ -62,12 +62,12 @@ def scalar_newton(fn, x0, atol = 1.0e-6, miter = 100):
         if torch.all(torch.abs(R) < atol):
             break
         
-        x -= R / J
+        x = x - R / J
         
         R, J = fn(x)
     else:
         warnings.warn("Scalar implicit solve did not succeed.  Results may be inaccurate...")
-
+    
     return x
 
 def scalar_bisection_newton(fn, a, b, atol = 1.0e-6, miter = 100, biter = 10):
diff --git a/test/test_flowrules.py b/test/test_flowrules.py
index e9a9d9b..65c54dd 100644
--- a/test/test_flowrules.py
+++ b/test/test_flowrules.py
@@ -18,7 +18,7 @@ def test_flow_rate(self):
             self.s,
             nbatch_dim=1,
         )
-        
+
         self.assertTrue(np.allclose(exact, numer, rtol=self.rtol))
 
     def test_history_rate(self):
@@ -46,7 +46,7 @@ def test_flow_history(self):
             self.h,
             nbatch_dim=1,
         ).unsqueeze(1)
-        
+
         self.assertTrue(np.allclose(exact, numer, rtol=self.rtol))
 
     def test_history_stress(self):
@@ -734,7 +734,7 @@ def setUp(self):
         self.kin = hardening.ChabocheHardeningModel(CP(self.C), CP(self.g))
 
         self.model = flowrules.IsoKinRateIndependentPlasticity(
-            CP(self.E), CP(self.sy), self.iso, self.kin
+            CP(self.E), CP(self.sy), self.iso, self.kin, s = 0.01
         )
 
         self.s = torch.linspace(160, 240, self.nbatch)
@@ -782,9 +782,33 @@ def test_d_residual_de(self):
         numer = utility.batch_differentiate(lambda x: self.model.ep_residual(self.ep_guess, self.s, self.h, self.t, self.T, x), self.erate).unsqueeze(-1)
         self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol))
 
-    def test_definiition_flow_rate(self):
-        vals = self.model.flow_rate(self.s, self.h, self.t, self.T, self.erate)[0]
-        yielding = self.model.f(self.s, self.h, self.T) > 0
-        r = self.model.ep_residual(vals, self.s, self.h, self.t, self.T, self.erate)
-        self.assertTrue(torch.allclose(r[yielding], torch.tensor(0.0)))
+    def test_d_sig(self):
+        x = torch.linspace(-1, 1, 10)
+        exact = self.model.dsig(x)
+        numer = utility.batch_differentiate(lambda x: self.model.sig(x), x)
+
+        self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol))
+
+    def test_df_ds(self):
+        exact = self.model.df_ds(self.s, self.h, self.T)
+        numer = utility.batch_differentiate(lambda x: self.model.f(x, self.h, self.T), self.s)
+
+        self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol))
+
+    def test_dmix_fn_ds(self):
+        exact = self.model.dmix_fn_ds(self.s, self.h, self.T)
+        numer = utility.batch_differentiate(lambda x: self.model.mix_fn(x, self.h, self.T), self.s)
+        
+        self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol, atol = 1.0e-4))
+
+    def test_df_dh(self):
+        exact = self.model.df_dh(self.s, self.h, self.T)
+        numer = utility.batch_differentiate(lambda x: self.model.f(self.s, x, self.T), self.h)
 
+        self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol))
+
+    def test_dmix_fn_dh(self):
+        exact = self.model.dmix_fn_dh(self.s, self.h, self.T)
+        numer = utility.batch_differentiate(lambda x: self.model.mix_fn(self.s, x, self.T), self.h)
+        
+        self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol))
diff --git a/test/test_models.py b/test/test_models.py
index 61e87a0..148a56d 100644
--- a/test/test_models.py
+++ b/test/test_models.py
@@ -71,7 +71,12 @@ def test_derivs_stress(self):
         ddv = utility.batch_differentiate(
             lambda x: use.forward(self.t, x)[0], self.state_stress
         )
-
+        
+        # 3 time steps, 3 state variables
+        print(stress_rates.shape)
+        print(dv[1])
+        print(ddv[1])
+        
         self.assertTrue(np.allclose(dv, ddv, rtol=1e-4, atol=1e-4))
 
     def test_partial_state(self):
@@ -114,6 +119,9 @@ def test_partial_erate(self):
             lambda y: self.model.forward(t, self.state_strain, y, T)[0], erate
         ).squeeze(-1)
 
+        print(exact)
+        print(numer)
+
         self.assertTrue(torch.allclose(exact, numer, atol=1e-4, rtol=1e-4))
 
 
@@ -330,13 +338,13 @@ def setUp(self):
 class TestIsoKinRIPlasticity(unittest.TestCase, CommonModel, CommonModelBatchBatch):
     def setUp(self):
         self.E = torch.tensor(100000.0)
-        self.sy = torch.tensor(10.0)
+        self.sy = torch.tensor(100.0)
 
         self.R = torch.tensor(101.0)
         self.d = torch.tensor(1.3)
         self.iso = hardening.VoceIsotropicHardeningModel(CP(self.R), CP(self.d))
 
-        self.C = torch.tensor(12000.0)
+        self.C = torch.tensor(1200.0)
         self.g = torch.tensor(10.1)
         self.kin = hardening.FAKinematicHardeningModel(CP(self.C), CP(self.g))
 
@@ -364,7 +372,7 @@ def setUp(self):
         ) / 3.0
         self.state_stress = (
             torch.tensor([[0.05, 30.0, 10.0], [0.07, 10.0, 15.0], [0.08, 50.0, 60.0]])
-        ) / 3.0
+        )
 
         self.t = self.times[2]
         self.step = 2