diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 0f2475523..fd43add3b 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -12,7 +12,7 @@ variables:
     value: main
   SPACK_BRANCH:
     description: Branch of BlueBrain Spack to use for the CI pipeline
-    value: develop
+    value: jelic/nmodl_sympy
   BLUECONFIGS_BRANCH:
     description: Branch of blueconfigs to trigger the simulation stack pipeline from
     value: main
@@ -46,6 +46,8 @@ trigger cvf:
   variables:
     # Tell CVF to use the same commits/branches as NMODL.
     SPACK_ENV_FILE_URL: $SPACK_SETUP_COMMIT_MAPPING_URL
+    SPACK_BRANCH:
+      value: jelic/nmodl_sympy
 
 simulation_stack:
   stage: .pre
@@ -63,11 +65,15 @@ simulation_stack:
     # the blueconfigs CI as well to let it know about both CVF
     # and BLUECONFIGS
     SPACK_SETUP_IGNORE_PACKAGE_VARIABLES: "CVF BLUECONFIGS"
+    SPACK_BRANCH:
+      value: jelic/nmodl_sympy
 
 .spack_nmodl:
   variables:
     SPACK_PACKAGE: nmodl
     SPACK_PACKAGE_SPEC: +python+tests
+    SPACK_BRANCH:
+      value: jelic/nmodl_sympy
 
 spack_setup:
   extends: .spack_setup_ccache
@@ -85,6 +91,8 @@ spack_setup:
     PARSE_GITHUB_PR_DESCRIPTIONS: "true"
     # Ignore CVF ang BLUECONFIGS branches since those don't have a spack package
     SPACK_SETUP_IGNORE_PACKAGE_VARIABLES: "CVF BLUECONFIGS"
+    SPACK_BRANCH:
+      value: jelic/nmodl_sympy
 
 build:intel:
   extends:
@@ -92,6 +100,8 @@ build:intel:
     - .spack_nmodl
   variables:
     SPACK_PACKAGE_COMPILER: oneapi
+    SPACK_BRANCH:
+      value: jelic/nmodl_sympy
 
 build:nvhpc:
   extends:
@@ -100,11 +110,15 @@ build:nvhpc:
   variables:
     SPACK_PACKAGE_COMPILER: nvhpc
     SPACK_PACKAGE_DEPENDENCIES: ^bison%gcc^flex%gcc^py-jinja2%gcc^py-sympy%gcc^py-pyyaml%gcc
+    SPACK_BRANCH:
+      value: jelic/nmodl_sympy
 
 .nmodl_tests:
   variables:
     # https://github.com/BlueBrain/nmodl/issues/737
     bb5_ntasks: 1
+    SPACK_BRANCH:
+      value: jelic/nmodl_sympy
 
 test:intel:
   extends:
diff --git a/pyproject.toml b/pyproject.toml
index 28082e014..979189ace 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -9,7 +9,7 @@ readme = {file = "README.rst", content-type = "text/x-rst"}
 dynamic = ["version"]
 dependencies = [
     "find_libpython",
-    "sympy>=1.3",
+    "sympy>=1.13",
     "importlib-metadata;python_version<'3.9'",
     "importlib-resources;python_version<'3.9'",
 ]
diff --git a/python/nmodl/ode.py b/python/nmodl/ode.py
index e5cc926d1..ff54a013e 100644
--- a/python/nmodl/ode.py
+++ b/python/nmodl/ode.py
@@ -308,14 +308,14 @@ def solve_lin_system(
                 order="canonical",
             )
             for var, expr in sub_exprs:
-                new_local_vars.append(sp.ccode(var))
+                new_local_vars.append(sp.ccode(var, strict=True))
                 code.append(
-                    f"{var} = {sp.ccode(expr.evalf(), user_functions=custom_fcts)}"
+                    f"{var} = {sp.ccode(expr.evalf(), user_functions=custom_fcts, strict=True)}"
                 )
             solution_vector = simplified_solution_vector[0]
         for var, expr in zip(state_vars, solution_vector):
             code.append(
-                f"{sp.ccode(var)} = {sp.ccode(expr.evalf(), contract=False, user_functions=custom_fcts)}"
+                f"{sp.ccode(var, strict=True)} = {sp.ccode(expr.evalf(), contract=False, user_functions=custom_fcts, strict=True)}"
             )
     else:
         # large linear system: construct and return matrix J, vector F such that
@@ -326,7 +326,7 @@ def solve_lin_system(
         vecFcode = []
         for i, expr in enumerate(vecF):
             vecFcode.append(
-                f"F[{i}] = {sp.ccode(expr.simplify().evalf(), user_functions=custom_fcts)}"
+                f"F[{i}] = {sp.ccode(expr.simplify().evalf(), user_functions=custom_fcts, strict=True)}"
             )
         # construct matrix J
         vecJcode = []
@@ -334,7 +334,7 @@ def solve_lin_system(
             # todo: fix indexing to be ascending order
             flat_index = matJ.rows * (i % matJ.rows) + (i // matJ.rows)
             vecJcode.append(
-                f"J[{flat_index}] = {sp.ccode(expr.simplify().evalf(), user_functions=custom_fcts)}"
+                f"J[{flat_index}] = {sp.ccode(expr.simplify().evalf(), user_functions=custom_fcts, strict=True)}"
             )
         # interweave
         code = _interweave_eqs(vecFcode, vecJcode)
@@ -375,7 +375,7 @@ def solve_non_lin_system(eq_strings, vars, constants, function_calls):
     vecFcode = []
     for i, eq in enumerate(eqs):
         vecFcode.append(
-            f"F[{i}] = {sp.ccode(eq.simplify().subs(X_vec_map).evalf(), user_functions=custom_fcts)}"
+            f"F[{i}] = {sp.ccode(eq.simplify().subs(X_vec_map).evalf(), user_functions=custom_fcts, strict=True)}"
         )
 
     vecJcode = []
@@ -385,6 +385,7 @@ def solve_non_lin_system(eq_strings, vars, constants, function_calls):
         rhs = sp.ccode(
             jacobian[i, j].simplify().subs(X_vec_map).evalf(),
             user_functions=custom_fcts,
+            strict=True,
         )
         vecJcode.append(f"J[{flat_index}] = {rhs}")
 
@@ -500,7 +501,7 @@ def integrate2c(diff_string, dt_var, vars, use_pade_approx=False):
     # 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(), strict=True)}"
 
 
 def forwards_euler2c(diff_string, dt_var, vars, function_calls):
@@ -528,7 +529,9 @@ def forwards_euler2c(diff_string, dt_var, vars, function_calls):
 
     custom_fcts = _get_custom_functions(function_calls)
     # return result as C code in NEURON format
-    return f"{sp.ccode(x)} = {sp.ccode(solution, user_functions=custom_fcts)}"
+    return (
+        f"{sp.ccode(x)} = {sp.ccode(solution, user_functions=custom_fcts, strict=True)}"
+    )
 
 
 def differentiate2c(expression, dependent_var, vars, prev_expressions=None):
@@ -612,9 +615,9 @@ def differentiate2c(expression, dependent_var, vars, prev_expressions=None):
     try:
         # if expression is equal to one of the supplied vars, replace with this var
         # can do a simple string comparison here since a var cannot be further simplified
-        diff_as_string = sp.ccode(diff)
+        diff_as_string = sp.ccode(diff, strict=True)
         for v in sympy_vars:
-            if diff_as_string == sp.ccode(sympy_vars[v]):
+            if diff_as_string == sp.ccode(sympy_vars[v], strict=True):
                 diff = sympy_vars[v]
 
         # or if equal to rhs of one of the supplied equations, replace with lhs
@@ -635,4 +638,4 @@ def differentiate2c(expression, dependent_var, vars, prev_expressions=None):
         pass
 
     # return result as C code in NEURON format
-    return sp.ccode(diff.evalf())
+    return sp.ccode(diff.evalf(), strict=True)
diff --git a/requirements.txt b/requirements.txt
index dbdba4c48..f032d00f1 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -5,7 +5,7 @@ Jinja2>=2.9.3
 PyYAML>=3.13
 # runtime dependencies
 find_libpython
-sympy>=1.3
+sympy>=1.13
 importlib-metadata;python_version<'3.9'
 importlib-resources;python_version<'3.9'
 # dependencies for test
diff --git a/test/unit/codegen/codegen_neuron_cpp_visitor.cpp b/test/unit/codegen/codegen_neuron_cpp_visitor.cpp
index 8268eb22d..4576e5e48 100644
--- a/test/unit/codegen/codegen_neuron_cpp_visitor.cpp
+++ b/test/unit/codegen/codegen_neuron_cpp_visitor.cpp
@@ -245,3 +245,23 @@ SCENARIO("ARTIFICIAL_CELL with `net_move`") {
         }
     }
 }
+
+SCENARIO("Codegen fails when analytic solution not available",
+         "[codegen][solver][sympy][derivative]") {
+    GIVEN("DERIVATIVE block with an equation which is not analytically solvable") {
+        std::string nmodl_text = R"(
+            BREAKPOINT {
+              SOLVE equation METHOD cnexp
+            }
+
+            DERIVATIVE equation {
+            var1' = sin(var1 * var1)
+            })";
+        THEN("Run Kinetic and Sympy Visitor and throw error") {
+            const auto& ast = NmodlDriver().parse_string(nmodl_text);
+            REQUIRE_THROWS_WITH(transpile(nmodl_text),
+                                ContainsSubstring(
+                                    "PRIME encountered during code generation, ODEs not solved?"));
+        }
+    }
+}