diff --git a/python/nmodl/ode.py b/python/nmodl/ode.py
index c1c907eae..3ae19cc19 100644
--- a/python/nmodl/ode.py
+++ b/python/nmodl/ode.py
@@ -539,10 +539,12 @@ def integrate2c(diff_string, dt_var, vars, use_pade_approx=False):
         if _a1 == 0 and _a2 == 0:
             solution = _a0
 
+    custom_fcts = {str(f.func): str(f.func) for f in solution.atoms(sp.Function)}
+
     # return result as C code in NEURON format:
     #   - in the lhs x_0 refers to the state var at time (t+dt)
     #   - in the rhs x_0 refers to the state var at time t
-    return f"{sp.ccode(x)} = {sp.ccode(solution.evalf())}"
+    return f"{sp.ccode(x)} = {sp.ccode(solution.evalf(), user_functions=custom_fcts)}"
 
 
 def forwards_euler2c(diff_string, dt_var, vars, function_calls):
diff --git a/test/unit/ode/test_ode.py b/test/unit/ode/test_ode.py
index 82e0358d2..6508fc0c2 100644
--- a/test/unit/ode/test_ode.py
+++ b/test/unit/ode/test_ode.py
@@ -169,6 +169,10 @@ def test_integrate2c():
         ("a", "x + a*dt"),
         ("a*x", "x*exp(a*dt)"),
         ("a*x+b", "(-b + (a*x + b)*exp(a*dt))/a"),
+        (
+            "a*x+custom_function(a)",
+            "(-custom_function(a) + (a*x + custom_function(a))*exp(a*dt))/a",
+        ),
     ]
     for eq, sol in test_cases:
         assert _equivalent(