diff --git a/tests/pools/common/unitary/test_remove_liquidity_one_coin.py b/tests/pools/common/unitary/test_remove_liquidity_one_coin.py
index 3a37d6ac..83dbfb6d 100644
--- a/tests/pools/common/unitary/test_remove_liquidity_one_coin.py
+++ b/tests/pools/common/unitary/test_remove_liquidity_one_coin.py
@@ -1,5 +1,6 @@
 import brownie
 import pytest
+from simulation import Curve
 
 pytestmark = [
     pytest.mark.skip_pool("busd", "compound", "pax", "susd", "usdt", "y"),
@@ -81,3 +82,13 @@ def test_below_zero(alice, swap):
 def test_above_n_coins(alice, swap, wrapped_coins, n_coins):
     with brownie.reverts():
         swap.remove_liquidity_one_coin(1, n_coins, 0, {'from': alice})
+
+@pytest.mark.itercoins("idx")
+def test_against_simulation(alice, swap, n_coins, pool_token, idx):
+    amount = pool_token.balanceOf(alice)
+
+    balances = [swap.balances(i) for i in range(n_coins)]
+    curve_model = Curve(swap.A(), balances, n_coins, tokens=pool_token.totalSupply())
+    expected, _ = curve_model.calc_withdraw_one_coin(amount, idx)
+
+    assert swap.calc_withdraw_one_coin(amount, idx) == expected
diff --git a/tests/simulation.py b/tests/simulation.py
index 5559b320..9fa1bddf 100644
--- a/tests/simulation.py
+++ b/tests/simulation.py
@@ -148,13 +148,24 @@ def remove_liquidity_imbalance(self, amounts):
 
     def calc_withdraw_one_coin(self, token_amount, i):
         xp = self.xp()
-        if self.fee:
-            fee = self.fee - self.fee * xp[i] // sum(xp) + 5 * 10 ** 5
-        else:
-            fee = 0
+        xp_reduced = list(xp)
 
         D0 = self.D()
         D1 = D0 - token_amount * D0 // self.tokens
-        dy = xp[i] - self.y_D(i, D1)
-
-        return dy - dy * fee // 10 ** 10
+        new_y = self.y_D(i, D1)
+
+        fee = self.fee * self.n // (4 * (self.n - 1))
+        for j in range(self.n):
+          dx_expected = 0
+          if j == i:
+            dx_expected = xp[j] * D1 // D0 - new_y
+          else:
+            dx_expected = xp[j] - xp[j] * D1 // D0
+          xp_reduced[j] -= fee * dx_expected // 10 ** 10
+
+        self.x = [x // (p // 10 ** 18) for x, p in zip(xp_reduced, self.p)]
+        dy = xp_reduced[i] - self.y_D(i, D1) - 1    # Withdraw less to account for rounding errors
+        self.x = [x // (p // 10 ** 18) for x, p in zip(xp, self.p)]
+        dy_0 = xp[i] - new_y
+
+        return dy, dy_0 - dy