diff --git a/src/nmreval/fit/_meta.py b/src/nmreval/fit/_meta.py
index 9fe9c34..889cc57 100644
--- a/src/nmreval/fit/_meta.py
+++ b/src/nmreval/fit/_meta.py
@@ -24,9 +24,10 @@ class ModelFactory:
             param_len.append(len(func['func'].params))
 
             if func['children']:
-                right, _, _ = ModelFactory.create_from_list(func['children'], left=func['func'], left_cnt=func['pos'],
+                right, _, _ = ModelFactory.create_from_list(func['children'], left_cnt=func['pos'],
                                                             func_order=func_order, param_len=param_len)
                 right_cnt = None
+                right = MultiModel(func['func'], right, func['children'][0]['op'], left_idx=func['cnt'], right_idx=None)
             else:
                 right = func['func']
                 right_cnt = func['cnt']
@@ -46,7 +47,13 @@ class MultiModel:
     str_op = {'+': operator.add, '*': operator.mul, '-': operator.sub, '/': operator.truediv}
     int_op = {0: operator.add, 1: operator.mul, 2: operator.sub, 3: operator.truediv}
 
-    def __init__(self, left: Any, right: Any, op: str | Callable | int = '+', left_idx=0, right_idx=1):
+    def __init__(self,
+                 left: Any,
+                 right: Any,
+                 op: str | Callable | int = '+',
+                 left_idx: int | None = 0,
+                 right_idx: int | None = 1,
+                 ):
         self._left = left
         self._right = right
 
diff --git a/src/nmreval/models/diffusion.py b/src/nmreval/models/diffusion.py
index 4e03f92..f99e97a 100644
--- a/src/nmreval/models/diffusion.py
+++ b/src/nmreval/models/diffusion.py
@@ -125,7 +125,7 @@ class Peschier:
         q = nucleus*g*tp
 
         r1s, r1f = 1 / t1s, 1 / t1f
-        kf, ks = pf*k, (1-pf)*k
+        kf, ks = (1-pf)*k, pf*k
         a_plus = 0.5 * (d*q*q + kf + ks + r1f + r1s + np.sqrt((d*q*q + kf + r1f - ks - r1s)**2 + 4*kf*ks))
         a_minu = 0.5 * (d*q*q + kf + ks + r1f + r1s - np.sqrt((d*q*q + kf + r1f - ks - r1s)**2 + 4*kf*ks))
 
@@ -149,7 +149,6 @@ class DiffusionGradients:
         else:
             tm = t2
             tp = x
-            # T2 decay happens twice
         q2 = (g1**2 - g2**2) * (nucleus * tp)**2
         t_eff = (2 * tp) / 3 + tm