diff --git a/phasefield/phasefieldstep.py b/phasefield/phasefieldstep.py
index 26d6121d20f5f926d43d9dd6facfda8932c5258e..706e0c0776f581bcb390b072d85c4cddc3e01ca7 100644
--- a/phasefield/phasefieldstep.py
+++ b/phasefield/phasefieldstep.py
@@ -121,7 +121,7 @@ class PhaseFieldStep:
         extra_force = sp.Matrix([0] * self.data_handling.dim)
         if order_parameter_force is not None:
             for order_parameter_idx, force in order_parameter_force.items():
-                extra_force += self.phi_field(order_parameter_idx) * sp.Matrix(force)
+                extra_force += sp.Matrix([f_i * self.phi_field(order_parameter_idx) for f_i in force])
         self.force_eqs = force_kernel_using_pressure_tensor(self.force_field, self.pressure_tensor_field, dx=dx,
                                                             extra_force=extra_force, discretization=discretization)
         self.force_from_pressure_tensor_kernel = create_kernel(apply_neumann_boundaries(self.force_eqs),