diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc
index 75ae03e88001d0d68d17a6faf1b0a973e378e18e..aebdf3210b3e1b3fd33d43002cf255eb20c7ccc1 100644
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2003-2020 Dynare Team
+ * Copyright © 2003-2021 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -543,75 +543,6 @@ ForecastStatement::writeJsonOutput(ostream &output) const
   output << "}";
 }
 
-DetCondForecastStatement::DetCondForecastStatement(SymbolList symbol_list_arg,
-                                                   OptionsList options_list_arg,
-                                                   const bool linear_decomposition_arg) :
-  options_list{move(options_list_arg)},
-  symbol_list{move(symbol_list_arg)},
-  linear_decomposition{linear_decomposition_arg}
-{
-
-}
-
-void
-DetCondForecastStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
-{
-  options_list.writeOutput(output);
-  if (linear_decomposition)
-    {
-      output << "first_order_solution_to_compute = 1;" << endl
-             << "if isfield(oo_, 'dr')" << endl
-             << "    if isfield(oo_.dr, 'ghx') && isfield(oo_.dr, 'ghu') && isfield(oo_.dr, 'state_var') && isfield(oo_.dr, 'order_var')" << endl
-             << "        first_order_solution_to_compute = 0;" << endl
-             << "    end" << endl
-             << "end" << endl
-             << "if first_order_solution_to_compute" << endl
-             << "  fprintf('%s','Computing the first order solution ...');" << endl
-             << "  options_.nograph = true;" << endl
-             << "  options_.order = 1;" << endl
-             << "  options_.noprint = true;" << endl
-             << "  options_.nocorr = true;" << endl
-             << "  options_.nomoments = true;" << endl
-             << "  options_.nodecomposition = true;" << endl
-             << "  options_.nofunctions = true;" << endl
-             << "  options_.irf = 0;" << endl
-             << "  tmp_periods = options_.periods;" << endl
-             << "  options_.periods = 0;" << endl
-             << "  var_list_ = char();" << endl
-             << "  info = stoch_simul(var_list_);" << endl
-             << R"(  fprintf('%s\n','done');)" << endl
-             << "  options_.periods = tmp_periods;" << endl
-             << "end" << endl;
-    }
-  vector<string> symbols = symbol_list.get_symbols();
-  if (symbols.size() > 0)
-    output << symbols[1] << " = det_cond_forecast(";
-  for (unsigned int i = 0; i < symbols.size() - 1; i++)
-    output << symbols[i] << ", ";
-  if (symbols.size() > 0)
-    output << symbols[symbols.size() - 1];
-  output << ");" << endl;
-}
-
-void
-DetCondForecastStatement::writeJsonOutput(ostream &output) const
-{
-  output << R"({"statementName": "det_cond_forecast")";
-  if (options_list.getNumberOfOptions())
-    {
-      output << ", ";
-      options_list.writeJsonOutput(output);
-    }
-  if (!symbol_list.empty())
-    {
-      output << ", ";
-      symbol_list.writeJsonOutput(output);
-    }
-  output << R"(, "linear_decomposition": )"
-         << (linear_decomposition ? "true" : "false")
-         << "}";
-}
-
 RamseyModelStatement::RamseyModelStatement(OptionsList options_list_arg) :
   options_list{move(options_list_arg)}
 {
diff --git a/src/ComputingTasks.hh b/src/ComputingTasks.hh
index 6bc0d05072b73f390c43d011a189f1d1dd5b75f2..aac12adec984f9e0d44882961c2a68ccccf78fac 100644
--- a/src/ComputingTasks.hh
+++ b/src/ComputingTasks.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2003-2020 Dynare Team
+ * Copyright © 2003-2021 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -95,20 +95,6 @@ public:
   void writeJsonOutput(ostream &output) const override;
 };
 
-class DetCondForecastStatement : public Statement
-{
-private:
-  const OptionsList options_list;
-  const SymbolList symbol_list;
-  const bool linear_decomposition;
-public:
-  DetCondForecastStatement(SymbolList symbol_list_arg,
-                           OptionsList options_list_arg,
-                           const bool linear_decompositiontion_arg);
-  void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
-  void writeJsonOutput(ostream &output) const override;
-};
-
 class ModelInfoStatement : public Statement
 {
 private:
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index c4573c5acd4bafdf9e4680ae6ffd3017ed0fd6f2..3582a97a1cfe368b793bb624dc39e923b133f069 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2003-2020 Dynare Team
+ * Copyright © 2003-2021 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -1032,7 +1032,7 @@ DynamicModel::writeDynamicBytecode(const string &basename) const
 }
 
 void
-DynamicModel::writeDynamicBlockBytecode(const string &basename, bool linear_decomposition) const
+DynamicModel::writeDynamicBlockBytecode(const string &basename) const
 {
   struct Uff_l
   {
@@ -1056,11 +1056,7 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename, bool linear_deco
   map<expr_t, int> reference_count;
   vector<int> feedback_variables;
   bool file_open = false;
-  string main_name;
-  if (linear_decomposition)
-    main_name = basename + "/model/bytecode/non_linear.cod";
-  else
-    main_name = basename + "/model/bytecode/dynamic.cod";
+  string main_name = basename + "/model/bytecode/dynamic.cod";
   code_file.open(main_name, ios::out | ios::binary | ios::ate);
   if (!code_file.is_open())
     {
@@ -1095,7 +1091,7 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename, bool linear_deco
           || simulation_type == BlockSimulationType::solveForwardComplete)
         {
           writeBlockBytecodeBinFile(basename, block, u_count_int, file_open,
-                                    simulation_type == BlockSimulationType::solveTwoBoundariesComplete || simulation_type == BlockSimulationType::solveTwoBoundariesSimple, linear_decomposition);
+                                    simulation_type == BlockSimulationType::solveTwoBoundariesComplete || simulation_type == BlockSimulationType::solveTwoBoundariesSimple);
           file_open = true;
         }
 
@@ -1123,8 +1119,6 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename, bool linear_deco
       fbeginblock.write(code_file, instruction_number);
 
       temporary_terms_t temporary_terms_union;
-      if (linear_decomposition)
-        compileTemporaryTerms(code_file, instruction_number, true, false, temporary_terms_union, blocks_temporary_terms_idxs);
 
       //The Temporary terms
       deriv_node_temp_terms_t tef_terms;
@@ -1155,8 +1149,7 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename, bool linear_deco
       // The equations
       for (i = 0; i < block_size; i++)
         {
-          if (!linear_decomposition)
-            write_eq_tt(i);
+          write_eq_tt(i);
 
           int variable_ID, equation_ID;
           EquationType equ_type;
@@ -1228,8 +1221,7 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename, bool linear_deco
           && simulation_type != BlockSimulationType::evaluateForward)
         {
           // Write temporary terms for derivatives
-          if (!linear_decomposition)
-            write_eq_tt(blocks[block].size);
+          write_eq_tt(blocks[block].size);
 
           switch (simulation_type)
             {
@@ -1556,16 +1548,11 @@ DynamicModel::printNonZeroHessianEquations(ostream &output) const
 
 void
 DynamicModel::writeBlockBytecodeBinFile(const string &basename, int num, int &u_count_int,
-                                        bool &file_open, bool is_two_boundaries, bool linear_decomposition) const
+                                        bool &file_open, bool is_two_boundaries) const
 {
   int j;
   std::ofstream SaveCode;
-  string filename;
-
-  if (!linear_decomposition)
-    filename = basename + "/model/bytecode/dynamic.bin";
-  else
-    filename = basename + "/model/bytecode/non_linear.bin";
+  string filename = basename + "/model/bytecode/dynamic.bin";
 
   if (file_open)
     SaveCode.open(filename, ios::out | ios::in | ios::binary | ios::ate);
@@ -2811,7 +2798,7 @@ DynamicModel::writeBlockDriverOutput(ostream &output, const string &basename, co
 }
 
 void
-DynamicModel::writeDriverOutput(ostream &output, const string &basename, bool block_decomposition, bool linear_decomposition, bool use_dll, bool estimation_present, bool compute_xrefs, bool julia) const
+DynamicModel::writeDriverOutput(ostream &output, const string &basename, bool block_decomposition, bool use_dll, bool estimation_present, bool compute_xrefs, bool julia) const
 {
   /* Writing initialisation for M_.lead_lag_incidence matrix
      M_.lead_lag_incidence is a matrix with as many columns as there are
@@ -2983,7 +2970,7 @@ DynamicModel::writeDriverOutput(ostream &output, const string &basename, bool bl
         }
 
   // Write the block structure of the model
-  if (block_decomposition || linear_decomposition)
+  if (block_decomposition)
     writeBlockDriverOutput(output, basename, modstruct, state_var, estimation_present);
 
   output << modstruct << "state_var = [";
@@ -4254,7 +4241,7 @@ DynamicModel::substitutePacExpectation(const string &pac_model_name)
 void
 DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsOrder,
                             const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll,
-                            bool bytecode, bool linear_decomposition)
+                            bool bytecode)
 {
   assert(jacobianExo || (derivsOrder < 2 && paramsDerivsOrder == 0));
 
@@ -4276,8 +4263,7 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
     }
 
   // Launch computations
-  cout << "Computing " << (linear_decomposition ? "nonlinear " : "")
-       << "dynamic model derivatives (order " << derivsOrder << ")." << endl;
+  cout << "Computing dynamic model derivatives (order " << derivsOrder << ")." << endl;
 
   computeDerivatives(derivsOrder, vars);
 
@@ -4292,32 +4278,7 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
     }
 
 
-  if (linear_decomposition)
-    {
-      auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
-      equationLinear(first_order_endo_derivatives);
-
-      auto contemporaneous_jacobian = evaluateAndReduceJacobian(eval_context);
-
-      if (!computeNaturalNormalization())
-        computeNonSingularNormalization(contemporaneous_jacobian);
-
-      select_non_linear_equations_and_variables();
-
-      equationTypeDetermination(first_order_endo_derivatives, 0);
-
-      reduceBlockDecomposition();
-
-      computeChainRuleJacobian();
-
-      determineLinearBlocks();
-
-      computeBlockDynJacobianCols();
-
-      if (!no_tmp_terms)
-        computeBlockTemporaryTerms();
-    }
-  else if (block)
+  if (block)
     {
       auto contemporaneous_jacobian = evaluateAndReduceJacobian(eval_context);
 
@@ -4614,7 +4575,7 @@ DynamicModel::computeBlockDynJacobianCols()
 }
 
 void
-DynamicModel::writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const
+DynamicModel::writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const
 {
   filesystem::path model_dir{basename};
   model_dir /= "model";
@@ -4623,20 +4584,10 @@ DynamicModel::writeDynamicFile(const string &basename, bool block, bool linear_d
   if (bytecode)
     filesystem::create_directories(model_dir / "bytecode");
 
-  if (linear_decomposition)
-    {
-      if (bytecode)
-        writeDynamicBlockBytecode(basename, linear_decomposition);
-      else
-        {
-          cerr << "'linear_decomposition' option requires the 'bytecode' option" << endl;
-          exit(EXIT_FAILURE);
-        }
-    }
-  else if (block)
+  if (block)
     {
       if (bytecode)
-        writeDynamicBlockBytecode(basename, linear_decomposition);
+        writeDynamicBlockBytecode(basename);
       else if (use_dll)
         {
           writeDynamicPerBlockCFiles(basename);
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 7752964f71365a63333d3330c64aeb9d1e85cea7..1c68bb874dcf07e202b16473f7f154dc066a1a8f 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2003-2020 Dynare Team
+ * Copyright © 2003-2021 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -145,12 +145,12 @@ private:
   //! Writes the per-block dynamic files of block decomposed model (C version)
   void writeDynamicPerBlockCFiles(const string &basename) const;
   //! Writes the code of the block-decomposed model in virtual machine bytecode
-  void writeDynamicBlockBytecode(const string &basename, bool linear_decomposition) const;
+  void writeDynamicBlockBytecode(const string &basename) const;
   //! Writes the code of the model in virtual machine bytecode
   void writeDynamicBytecode(const string &basename) const;
   //! Adds per-block information for bytecode simulation in a separate .bin file
   void writeBlockBytecodeBinFile(const string &basename, int num, int &u_count_int, bool &file_open,
-                                 bool is_two_boundaries, bool linear_decomposition) const;
+                                 bool is_two_boundaries) const;
 
   void writeSetAuxiliaryVariables(const string &basename, bool julia) const;
   void writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const;
@@ -298,9 +298,9 @@ public:
     \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files
   */
   void computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsOrder,
-                     const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode, bool linear_decomposition);
+                     const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode);
   //! Writes information about the dynamic model to the driver file
-  void writeDriverOutput(ostream &output, const string &basename, bool block, bool linear_decomposition, bool use_dll, bool estimation_present, bool compute_xrefs, bool julia) const;
+  void writeDriverOutput(ostream &output, const string &basename, bool block, bool use_dll, bool estimation_present, bool compute_xrefs, bool julia) const;
 
   //! Write JSON AST
   void writeJsonAST(ostream &output) const;
@@ -385,7 +385,7 @@ public:
   void substitutePacExpectation(const string &pac_model_name);
 
   //! Writes dynamic model file
-  void writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const;
+  void writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const;
   //! Writes file containing parameters derivatives
   void writeParamsDerivativesFile(const string &basename, bool julia) const;
 
diff --git a/src/DynareBison.yy b/src/DynareBison.yy
index e19e73255ba3a3a8f9ca7c5df664fd0653a881bf..44d14904c5f36f066561d0fcc8f6e1977e49f25c 100644
--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
@@ -1,6 +1,6 @@
 // -*- C++ -*-
 /*
- * Copyright © 2003-2020 Dynare Team
+ * Copyright © 2003-2021 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -80,7 +80,7 @@ class ParsingDriver;
 %token BVAR_REPLIC BYTECODE ALL_VALUES_REQUIRED PROPOSAL_DISTRIBUTION REALTIME VINTAGE
 %token CALIB_SMOOTHER CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR CUTOFF CYCLE_REDUCTION LOGARITHMIC_REDUCTION
 %token COMMA CONSIDER_ALL_ENDOGENOUS CONSIDER_ONLY_OBSERVED INITIAL_CONDITION_DECOMPOSITION
-%token DATAFILE FILE SERIES DET_COND_FORECAST DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
+%token DATAFILE FILE SERIES DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
 %token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT EXTENDED_PATH ENDOGENOUS_PRIOR EXPRESSION
 %token FILENAME DIRNAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS FIRST_SIMULATION_PERIOD LAST_OBS 
 %token SET_TIME OSR_PARAMS_BOUNDS KEEP_KALMAN_ALGO_IF_SINGULARITY_IS_DETECTED
@@ -313,7 +313,6 @@ statement : parameters
           | method_of_moments
           | shock_groups
           | init2shocks
-          | det_cond_forecast
           | var_expectation_model
           | compilation_setup
           | matched_moments
@@ -865,7 +864,6 @@ model_options : BLOCK { driver.block(); }
               | DIFFERENTIATE_FORWARD_VARS EQUAL '(' symbol_list ')' { driver.differentiate_forward_vars_some(); }
               | o_linear
               | PARALLEL_LOCAL_FILES EQUAL '(' parallel_local_filename_list ')'
-              | LINEAR_DECOMPOSITION { driver.linear_decomposition(); }
               | BALANCED_GROWTH_TEST_TOL EQUAL non_negative_number { driver.balanced_growth_test_tol($3); }
               ;
 
@@ -1352,16 +1350,6 @@ prior_posterior_function_options : o_function
                                  | o_sampling_draws
                                  ;
 
-det_cond_forecast : DET_COND_FORECAST '(' symbol ')' ';'
-                    { driver.det_cond_forecast_linear_decomposition($3); }
-                  | DET_COND_FORECAST '(' symbol COMMA symbol ')' ';'
-                    { driver.det_cond_forecast_linear_decomposition($3,$5); }
-                  | DET_COND_FORECAST '(' symbol COMMA simul_options_list ')' ';'
-                    { driver.det_cond_forecast_linear_decomposition($3); }
-                  | DET_COND_FORECAST '(' symbol COMMA symbol COMMA simul_options_list ')' ';'
-                    { driver.det_cond_forecast_linear_decomposition($3,$5); }
-                  ;
-
 simul : SIMUL ';'
         { driver.simul(); }
       | SIMUL '(' simul_options_list ')' ';'
diff --git a/src/DynareFlex.ll b/src/DynareFlex.ll
index 124f8c9625d677ce1f8926ea16ff61652a76189c..390e87bfd042dcbf9fcb9d4d6effa74a51963946 100644
--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
@@ -1,6 +1,6 @@
 /* -*- C++ -*- */
 /*
- * Copyright © 2003-2020 Dynare Team
+ * Copyright © 2003-2021 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -185,7 +185,6 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>smoother2histval {BEGIN DYNARE_STATEMENT; return token::SMOOTHER2HISTVAL;}
 <INITIAL>perfect_foresight_setup {BEGIN DYNARE_STATEMENT; return token::PERFECT_FORESIGHT_SETUP;}
 <INITIAL>perfect_foresight_solver {BEGIN DYNARE_STATEMENT; return token::PERFECT_FORESIGHT_SOLVER;}
-<INITIAL>det_cond_forecast {BEGIN DYNARE_STATEMENT; return token::DET_COND_FORECAST;}
 <INITIAL>compilation_setup {BEGIN DYNARE_STATEMENT; return token::COMPILATION_SETUP;}
 
 <DYNARE_STATEMENT>; {
@@ -813,7 +812,6 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK>use_dll {return token::USE_DLL;}
 <DYNARE_BLOCK>block {return token::BLOCK;}
 <DYNARE_BLOCK>bytecode {return token::BYTECODE;}
-<DYNARE_BLOCK>linear_decomposition {return token::LINEAR_DECOMPOSITION;}
 <DYNARE_BLOCK>all_values_required {return token::ALL_VALUES_REQUIRED;}
 <DYNARE_BLOCK>no_static {return token::NO_STATIC;}
 <DYNARE_BLOCK>differentiate_forward_vars {return token::DIFFERENTIATE_FORWARD_VARS;}
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 48b4489d32b40cf138993f07c864677a38350479..4b5814bf16fbc3f54b05295ace3e8e3b76e8e1a3 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -45,8 +45,6 @@ ModFile::ModFile(WarningConsolidation &warnings_arg)
                                        trend_component_model_table, var_model_table},
     orig_ramsey_dynamic_model{symbol_table, num_constants, external_functions_table,
                               trend_component_model_table, var_model_table},
-    non_linear_equations_dynamic_model{symbol_table, num_constants, external_functions_table,
-                                       trend_component_model_table, var_model_table},
     epilogue{symbol_table, num_constants, external_functions_table,
              trend_component_model_table, var_model_table},
     static_model{symbol_table, num_constants, external_functions_table},
@@ -193,17 +191,6 @@ ModFile::checkPass(bool nostrict, bool stochastic)
       cerr << "ERROR: In 'model' block, 'use_dll' option is not compatible with 'bytecode'" << endl;
       exit(EXIT_FAILURE);
     }
-  if (block && linear_decomposition)
-    {
-      cerr << "ERROR: In 'model' block, 'block' option is not compatible with 'linear_decomposition'" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  if (!bytecode && linear_decomposition)
-    {
-      cerr << "ERROR: For the moment in 'model' block, 'linear_decomposition' option is compatible only with 'bytecode' option" << endl;
-      exit(EXIT_FAILURE);
-    }
 
   if ((stochastic_statement_present || mod_file_struct.check_present || mod_file_struct.steady_present) && no_static)
     {
@@ -731,12 +718,6 @@ ModFile::computingPass(bool no_tmp_terms, OutputType output, int params_derivs_o
 
       // Compute static model and its derivatives
       static_model = static_cast<StaticModel>(dynamic_model);
-      if (linear_decomposition)
-        {
-          non_linear_equations_dynamic_model = dynamic_model;
-          non_linear_equations_dynamic_model.set_cutoff_to_zero();
-          non_linear_equations_dynamic_model.computingPass(true, 1, 0, global_eval_context, no_tmp_terms, block, use_dll, bytecode, linear_decomposition);
-        }
       if (!no_static)
         {
           if (mod_file_struct.stoch_simul_present
@@ -771,7 +752,7 @@ ModFile::computingPass(bool no_tmp_terms, OutputType output, int params_derivs_o
                 derivsOrder = 2;
               else if  (output == OutputType::third)
                 derivsOrder = 3;
-              dynamic_model.computingPass(true, derivsOrder, 0, global_eval_context, no_tmp_terms, block, use_dll, bytecode, linear_decomposition);
+              dynamic_model.computingPass(true, derivsOrder, 0, global_eval_context, no_tmp_terms, block, use_dll, bytecode);
             }
           else
             {
@@ -800,13 +781,13 @@ ModFile::computingPass(bool no_tmp_terms, OutputType output, int params_derivs_o
                   || mod_file_struct.estimation_analytic_derivation
                   || (mod_file_struct.GMM_present && (mod_file_struct.analytic_standard_errors_present || mod_file_struct.analytic_jacobian_present)))
                 paramsDerivsOrder = params_derivs_order;
-              dynamic_model.computingPass(true, derivsOrder, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, bytecode, linear_decomposition);
+              dynamic_model.computingPass(true, derivsOrder, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, bytecode);
               if (linear && mod_file_struct.ramsey_model_present)
-                orig_ramsey_dynamic_model.computingPass(true, 2, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, bytecode, linear_decomposition);
+                orig_ramsey_dynamic_model.computingPass(true, 2, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, bytecode);
             }
         }
       else // No computing task requested, compute derivatives up to 2nd order by default
-        dynamic_model.computingPass(true, 2, 0, global_eval_context, no_tmp_terms, block, use_dll, bytecode, linear_decomposition);
+        dynamic_model.computingPass(true, 2, 0, global_eval_context, no_tmp_terms, block, use_dll, bytecode);
 
       /* Check that the model is linear.
          FIXME: this check always passes if derivsOrder = 1, i.e. for a perfect
@@ -842,7 +823,7 @@ ModFile::computingPass(bool no_tmp_terms, OutputType output, int params_derivs_o
   // Compute epilogue derivatives (but silence standard output)
   streambuf *oldcout = cout.rdbuf();
   cout.rdbuf(nullptr);
-  epilogue.computingPass(true, 2, 0, global_eval_context, true, false, false, false, false);
+  epilogue.computingPass(true, 2, 0, global_eval_context, true, false, false, false);
   cout.rdbuf(oldcout);
 }
 
@@ -989,8 +970,7 @@ ModFile::writeMOutput(const string &basename, bool clear_all, bool clear_global,
   mOutputFile << "options_.linear = " << to_matlab_logical(linear) << ";" << endl
               << "options_.block = " << to_matlab_logical(block) << ";" << endl
               << "options_.bytecode = " << to_matlab_logical(bytecode) << ";" << endl
-              << "options_.use_dll = " << to_matlab_logical(use_dll) << ";" << endl
-              << "options_.linear_decomposition = " << to_matlab_logical(linear_decomposition) << ";" << endl;
+              << "options_.use_dll = " << to_matlab_logical(use_dll) << ";" << endl;
 
   if (parallel_local_files.size() > 0)
     {
@@ -1032,9 +1012,7 @@ ModFile::writeMOutput(const string &basename, bool clear_all, bool clear_global,
 
   if (dynamic_model.equation_number() > 0)
     {
-      if (linear_decomposition)
-        non_linear_equations_dynamic_model.writeDriverOutput(mOutputFile, basename, block, true, use_dll, mod_file_struct.estimation_present, compute_xrefs, false);
-      dynamic_model.writeDriverOutput(mOutputFile, basename, block, false, use_dll, mod_file_struct.estimation_present, compute_xrefs, false);
+      dynamic_model.writeDriverOutput(mOutputFile, basename, block, use_dll, mod_file_struct.estimation_present, compute_xrefs, false);
       if (!no_static)
         static_model.writeDriverOutput(mOutputFile, block);
     }
@@ -1141,13 +1119,7 @@ ModFile::writeMOutput(const string &basename, bool clear_all, bool clear_global,
               static_model.writeParamsDerivativesFile(basename, false);
             }
 
-          if (linear_decomposition)
-            {
-              non_linear_equations_dynamic_model.writeDynamicFile(basename, block, linear_decomposition, bytecode, use_dll, mexext, matlabroot, dynareroot, false);
-              non_linear_equations_dynamic_model.writeParamsDerivativesFile(basename, false);
-            }
-
-          dynamic_model.writeDynamicFile(basename, block, false, bytecode, use_dll, mexext, matlabroot, dynareroot, false);
+          dynamic_model.writeDynamicFile(basename, block, bytecode, use_dll, mexext, matlabroot, dynareroot, false);
 
           dynamic_model.writeParamsDerivativesFile(basename, false);
 
@@ -1247,15 +1219,14 @@ ModFile::writeJuliaOutput(const string &basename) const
 
   if (dynamic_model.equation_number() > 0)
     {
-      dynamic_model.writeDriverOutput(jlOutputFile, basename, false, false, false,
+      dynamic_model.writeDriverOutput(jlOutputFile, basename, false, false,
                                       mod_file_struct.estimation_present, false, true);
       if (!no_static)
         {
           static_model.writeStaticFile(basename, false, false, false, "", {}, {}, true);
           static_model.writeParamsDerivativesFile(basename, true);
         }
-      dynamic_model.writeDynamicFile(basename, block, linear_decomposition, bytecode, use_dll,
-                                     "", {}, {}, true);
+      dynamic_model.writeDynamicFile(basename, block, bytecode, use_dll, "", {}, {}, true);
       dynamic_model.writeParamsDerivativesFile(basename, true);
     }
   steady_state_model.writeSteadyStateFile(basename, mod_file_struct.ramsey_model_present, true);
diff --git a/src/ModFile.hh b/src/ModFile.hh
index b2ba54c703769451d736f33eeaa593d1d7999b3c..01e1a471d1192d588533e5648efa8cdb7cb18740 100644
--- a/src/ModFile.hh
+++ b/src/ModFile.hh
@@ -67,8 +67,6 @@ public:
   DynamicModel ramsey_FOC_equations_dynamic_model;
   //! A copy of the original model, used to test model linearity under ramsey problem
   DynamicModel orig_ramsey_dynamic_model;
-  //! A copy of the Dynamic model, containing only the non-linear equations w.r. to endogenous variables
-  DynamicModel non_linear_equations_dynamic_model;
   //! Epilogue model, as declared in the "epilogue" block
   Epilogue epilogue;
   //! Static model, as derived from the "model" block when leads and lags have been removed
@@ -101,9 +99,6 @@ public:
     with a lead */
   vector<string> differentiate_forward_vars_subset;
 
-  //! Is the model is block-decomposed according the linear and the non-linear equations
-  bool linear_decomposition{false};
-
   //! Are nonstationary variables present ?
   bool nonstationary_variables{false};
 
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 6169fb70602c9b52a9d38e9e78d9e82c1af223fb..58def010de4f691e6fccb65b93cdbbcce4c0c212 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -351,35 +351,6 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context) const
   return contemporaneous_jacobian;
 }
 
-void
-ModelTree::select_non_linear_equations_and_variables()
-{
-  endo2block.resize(endo2eq.size(), 1); // The 1 is a dummy value, distinct from 0
-  eq2block.resize(endo2eq.size(), 1);
-  int i = 0;
-  for (int endo = 0; endo < static_cast<int>(endo2eq.size()); endo++)
-    {
-      int eq = endo2eq[endo];
-      if (!is_equation_linear[eq])
-        {
-          eq_idx_block2orig[i] = eq;
-          endo_idx_block2orig[i] = endo;
-          endo2block[endo] = 0;
-          eq2block[eq] = 0;
-          i++;
-        }
-    }
-  updateReverseVariableEquationOrderings();
-
-  blocks.clear();
-  blocks.resize(1);
-  blocks[0].size = i;
-  blocks[0].mfs_size = i;
-  blocks[0].first_equation = 0;
-  computeDynamicStructureOfBlock(0);
-  computeSimulationTypeOfBlock(0);
-}
-
 bool
 ModelTree::computeNaturalNormalization()
 {
@@ -870,23 +841,6 @@ ModelTree::reduceBlockDecomposition()
       }
 }
 
-void
-ModelTree::equationLinear(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives)
-{
-  is_equation_linear.clear();
-  is_equation_linear.resize(symbol_table.endo_nbr(), true);
-  for (const auto &[indices, expr] : first_order_endo_derivatives)
-    {
-      set<pair<int, int>> endogenous;
-      expr->collectEndogenous(endogenous);
-      if (endogenous.size() > 0)
-        {
-          int eq = get<0>(indices);
-          is_equation_linear[eq] = false;
-        }
-    }
-}
-
 void
 ModelTree::determineLinearBlocks()
 {
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 0ea39d5c46f5f499c13ec060a652588e5ae05eff..c65a8c3c068fdd8310040cf748a2d8c2e6e788fa 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2003-2020 Dynare Team
+ * Copyright © 2003-2021 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -302,8 +302,6 @@ protected:
   /*! Returns the contemporaneous_jacobian.
       Elements below the cutoff are discarded. External functions are evaluated to 1. */
   jacob_map_t evaluateAndReduceJacobian(const eval_context_t &eval_context) const;
-  //! Select and reorder the non linear equations of the model
-  void select_non_linear_equations_and_variables();
   /* Search the equations and variables belonging to the prologue and the
      epilogue of the model.
      Initializes “eq_idx_block2orig” and “endo_idx_block2orig”.
@@ -334,8 +332,6 @@ protected:
      max_lead) across all its occurences inside the equations of the block to
      which it belongs. */
   pair<lag_lead_vector_t, lag_lead_vector_t> getVariableLeadLagByBlock() const;
-  //! For each equation determine if it is linear or not
-  void equationLinear(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives);
   //! Print an abstract of the block structure of the model
   void printBlockDecomposition() const;
   //! Determine for each block if it is linear or not
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index aadae940d41ae37f41ae83250f7d48c565dc2c35..e95181174ea7c709da2b43b301280a721243b9e3 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -622,12 +622,6 @@ ParsingDriver::block()
   mod_file->block = true;
 }
 
-void
-ParsingDriver::linear_decomposition()
-{
-  mod_file->linear_decomposition = true;
-}
-
 void
 ParsingDriver::no_static()
 {
@@ -2348,27 +2342,6 @@ ParsingDriver::conditional_forecast_paths()
   det_shocks.clear();
 }
 
-void
-ParsingDriver::det_cond_forecast_linear_decomposition(const string &plan)
-{
-  symbol_list.clear();
-  symbol_list.addSymbol(plan);
-  mod_file->addStatement(make_unique<DetCondForecastStatement>(symbol_list, options_list, mod_file->linear_decomposition));
-  symbol_list.clear();
-  options_list.clear();
-}
-
-void
-ParsingDriver::det_cond_forecast_linear_decomposition(const string &plan, const string &dset)
-{
-  symbol_list.clear();
-  symbol_list.addSymbol(plan);
-  symbol_list.addSymbol(dset);
-  mod_file->addStatement(make_unique<DetCondForecastStatement>(symbol_list, options_list, mod_file->linear_decomposition));
-  symbol_list.clear();
-  options_list.clear();
-}
-
 void
 ParsingDriver::calib_smoother()
 {
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index 6f735acad1ee1c8954ecc76b0e84860b03eecdb8..04b791b4d10c2665102e984559e1f19a72f9aae7 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2003-2020 Dynare Team
+ * Copyright © 2003-2021 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -326,8 +326,6 @@ public:
   void use_dll();
   //! the modelis block decomposed
   void block();
-  //! the model is decomposed according to the linearity of its equations
-  void linear_decomposition();
 
   //! the model is stored in a binary file
   void bytecode();
@@ -697,9 +695,6 @@ public:
   void conditional_forecast_paths();
   //! Plot conditional forecast statement
   void plot_conditional_forecast(const string &periods = "");
-  //! Deterministic conditional forecast statement
-  void det_cond_forecast_linear_decomposition(const string &plan);
-  void det_cond_forecast_linear_decomposition(const string &plan, const string &dset);
   //! Smoother on calibrated models
   void calib_smoother();
   //! Extended path