diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp
index 4db3a9531..729b5865e 100644
--- a/src/math/lp/lar_solver.cpp
+++ b/src/math/lp/lar_solver.cpp
@@ -476,6 +476,9 @@ namespace lp {
         auto& x = m_mpq_lar_core_solver.m_r_x[j];
         auto delta = new_val - x;
         x = new_val;
+        TRACE("lar_solver_feas", tout << "setting " << j << " to "
+         << new_val << (column_is_feasible(j)?"feas":"non-feas") << "\n";);
+        m_mpq_lar_core_solver.m_r_solver.track_column_feasibility(j);
         change_basic_columns_dependend_on_a_given_nb_column(j, delta);
     }
 
@@ -1101,7 +1104,6 @@ namespace lp {
 
     mpq lar_solver::get_value(column_index const& j) const {
         SASSERT(get_status() == lp_status::OPTIMAL || get_status() == lp_status::FEASIBLE);
-        SASSERT(m_columns_with_changed_bounds.empty());
         numeric_pair<mpq> const& rp = get_column_value(j);
         return from_model_in_impq_to_mpq(rp);        
     }
@@ -1818,7 +1820,7 @@ namespace lp {
 
         if (is_base(j) && column_is_fixed(j))
             m_fixed_base_var_set.insert(j);
-        TRACE("lar_solver_feas", tout << "j = " << j << " became " << (this->column_is_feasible(j) ? "feas" : "non-feas") << ", and " << (this->column_is_bounded(j) ? "bounded" : "non-bounded") << std::endl;);    
+        TRACE("lar_solver_feas", tout << "j = " << j << " became " << (this->column_is_feasible(j) ? "feas" : "non-feas") << ", and " << (this->column_is_bounded(j) ? "bounded" : "non-bounded") << " val = " << get_column_value(j) << std::endl;);   
     }
 
     void lar_solver::insert_to_columns_with_changed_bounds(unsigned j) {
@@ -2363,6 +2365,11 @@ namespace lp {
         SASSERT(bdep != nullptr);
         m_crossed_bounds_deps = m_dependencies.mk_join(bdep, dep);
     }
+
+    void lar_solver::track_column_feasibility(lpvar j) {
+        m_mpq_lar_core_solver.m_r_solver.track_column_feasibility(j);
+    }
+
 } // namespace lp
 
 
diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h
index 2a1f8644c..0748d2507 100644
--- a/src/math/lp/lar_solver.h
+++ b/src/math/lp/lar_solver.h
@@ -92,6 +92,8 @@ class lar_solver : public column_namer {
 public: 
     const indexed_uint_set& columns_with_changed_bounds() const { return m_columns_with_changed_bounds; }
     inline void clear_columns_with_changed_bounds() { m_columns_with_changed_bounds.reset(); }
+    void track_column_feasibility(lpvar j);
+
 private:
     // m_touched_rows contains rows that have been changed by a pivoting or have a column with changed bounds
     indexed_uint_set m_touched_rows;
diff --git a/src/math/lp/lp_core_solver_base.h b/src/math/lp/lp_core_solver_base.h
index d7d1666d0..1d6634c06 100644
--- a/src/math/lp/lp_core_solver_base.h
+++ b/src/math/lp/lp_core_solver_base.h
@@ -571,7 +571,7 @@ public:
     }
     void remove_column_from_inf_heap(unsigned j) {
 		if (m_inf_heap.contains(j)) {
-            TRACE("lar_solver_inf_heap", tout << "insert into heap j = " << j << "\n";);
+            TRACE("lar_solver_inf_heap", tout << "erase from heap j = " << j << "\n";);
         	m_inf_heap.erase(j);
         }
         lp_assert(column_is_feasible(j));
diff --git a/src/math/lp/lp_core_solver_base_def.h b/src/math/lp/lp_core_solver_base_def.h
index 0552b8e82..b4a53e872 100644
--- a/src/math/lp/lp_core_solver_base_def.h
+++ b/src/math/lp/lp_core_solver_base_def.h
@@ -115,10 +115,12 @@ pretty_print(std::ostream & out) {
 
 template <typename T, typename X> void lp_core_solver_base<T, X>::
 add_delta_to_entering(unsigned entering, const X& delta) {
-    m_x[entering] += delta;     
+    m_x[entering] += delta;  
+    TRACE("lar_solver_feas", tout << "not tracking feas entering = " << entering << " = " << m_x[entering] << (column_is_feasible(entering) ? " feas" : " non-feas") << "\n";); 
     for (const auto & c : m_A.m_columns[entering]) {
         unsigned i = c.var();
         m_x[m_basis[i]] -= delta * m_A.get_val(c);
+        TRACE("lar_solver_feas", tout << "not tracking feas m_basis[i] = " << m_basis[i] << " = " << m_x[m_basis[i]] << (column_is_feasible(m_basis[i]) ? " feas" : " non-feas") << "\n";);
     }
 }
 
diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h
index 9871ec691..f5d7314db 100644
--- a/src/math/lp/lp_primal_core_solver.h
+++ b/src/math/lp/lp_primal_core_solver.h
@@ -394,9 +394,10 @@ namespace lp {
         const X &new_val_for_leaving = get_val_for_leaving(leaving);
         X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent;
         this->m_x[leaving] = new_val_for_leaving;
-        // this will remove the leaving from the heap
-        TRACE("lar_solver_inf_heap", tout << "leaving = " << leaving
+        TRACE("lar_solver_feas", tout << "entering = " << entering << ", leaving = " << leaving << ", new_val_for_leaving = " << new_val_for_leaving << ", theta = " << theta << "\n";);
+        TRACE("lar_solver_feas", tout << "leaving = " << leaving
                                  << " removed from inf_heap()\n";);
+        // this will remove the leaving from the heap
         this->inf_heap().erase_min();
         advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta);
         if (this->current_x_is_feasible())
diff --git a/src/math/lp/monomial_bounds.cpp b/src/math/lp/monomial_bounds.cpp
index 981677ab0..177cfb3d3 100644
--- a/src/math/lp/monomial_bounds.cpp
+++ b/src/math/lp/monomial_bounds.cpp
@@ -21,8 +21,14 @@ namespace nla {
         coeffs.push_back(std::make_pair(rational::one(), monic_var));
         lp::lpvar term_index = c().lra.add_term(coeffs, UINT_MAX);
         auto* dep = explain_fixed(vars, non_fixed);
+        // term_index becomes the column index of the term slack variable
         term_index = c().lra.map_term_index_to_column_index(term_index);
         c().lra.update_column_type_and_bound(term_index, lp::lconstraint_kind::EQ, mpq(0), dep);
+        c().lra.track_column_feasibility(term_index);
+        if (!c().lra.column_is_feasible(term_index)) {
+            c().lra.set_status(lp::lp_status::UNKNOWN);
+        }
+        
     }
 
     u_dependency* monomial_bounds::explain_fixed(const svector<lpvar>& vars, lpvar non_fixed) {
diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp
index 6641347e4..5a6af5659 100644
--- a/src/math/lp/nla_core.cpp
+++ b/src/math/lp/nla_core.cpp
@@ -138,6 +138,7 @@ void core::add_monic(lpvar v, unsigned sz, lpvar const* vs) {
         m_add_buffer[i] = j;
     }
     m_emons.add(v, m_add_buffer);
+    m_monics_with_changed_bounds.insert(v);
 }
     
 void core::push() {
@@ -1938,6 +1939,7 @@ void core::add_lower_bound_monic(lpvar j, const lp::mpq& v, bool is_strict, std:
     }
 
     void core::calculate_implied_bounds_for_monic(lp::lpvar monic_var) {
+        if (!is_monic_var(monic_var)) return;
         m_propagated.reserve(monic_var + 1, false);
         bool throttle = params().arith_nl_throttle_unit_prop();
         if (throttle && m_propagated[monic_var])
diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp
index d355b9118..2a7412c29 100644
--- a/src/smt/theory_lra.cpp
+++ b/src/smt/theory_lra.cpp
@@ -3191,6 +3191,8 @@ public:
     lbool make_feasible() {
         TRACE("pcs",  tout << lp().constraints(););
         TRACE("arith_verbose", tout << "before calling lp().find_feasible_solution()\n"; display(tout););
+        // todo: remove later : debug!!!!!
+        lp().move_non_basic_columns_to_bounds(false);
         auto status = lp().find_feasible_solution();
         TRACE("arith_verbose", display(tout););
         if (lp().is_feasible())