diff --git a/classical/systems.py b/classical/systems.py
index 230f6520f27631f46375570ea54f93efa64c2c31..6aae14796f058cf9d75a78980177512b3f951577 100644
--- a/classical/systems.py
+++ b/classical/systems.py
@@ -184,11 +184,11 @@ class ManyParticlePhysicalSystem:
 
     def get_potential_energy(self, qval):
         """ evaluate the potential energy numerically """
-        return float(self.potential_energy.subs(list(zip(self.q, qval))))
+        return float(self.potential_energy.xreplace(dict(zip(self.q, qval))))
 
     def get_kinetic_energy(self, pval):
         """ evaluate the kinetic energy numerically """
-        return float(self.kinetic_energy.subs(list(zip(self.p, pval))))
+        return float(self.kinetic_energy.xreplace(dict(zip(self.p, pval))))
 
     def get_total_energy(self, qval, pval):
         """ evaluate the total energy numerically """
@@ -196,7 +196,7 @@ class ManyParticlePhysicalSystem:
 
     def get_forces_subs(self, q):
         """ evaluate the forces numerically, subs version, very slow """
-        return self.forces.subs(list(zip(self.q, q)))
+        return self.forces.xreplace(dict(zip(self.q, q)))
 
     def get_forces_nopbc(self, qval):
         return np.array(self.get_forces_val(qval), dtype=float)
@@ -223,7 +223,7 @@ class MPPS1D(ManyParticlePhysicalSystem):
     def find_mic(self, qval):
         """ works for only two particles """
         micval = qval.copy()
-        r = float(self.ff[0].r_cart.subs(list(zip(self.q, qval))))
+        r = float(self.ff[0].r_cart.xreplace(dict(zip(self.q, qval))))
         r -= np.rint(r*self.inv_cell)*self.cell
         micval[1] = qval[0] + np.abs(r)
         return micval, np.abs(r)