diff --git a/russell_ode/src/bin/brusselator_pde.rs b/russell_ode/src/bin/brusselator_pde.rs
index c3958a4e..6aa15b37 100644
--- a/russell_ode/src/bin/brusselator_pde.rs
+++ b/russell_ode/src/bin/brusselator_pde.rs
@@ -1,6 +1,6 @@
 use russell_lab::{format_scientific, get_num_threads, set_num_threads, StrError};
 use russell_ode::prelude::*;
-use russell_sparse::Genie;
+use russell_sparse::{Genie, LinSolParams, Ordering, Scaling};
 use structopt::StructOpt;
 
 /// Command line options
@@ -38,6 +38,18 @@ struct Options {
     /// Use MUMPS solver
     #[structopt(long)]
     mumps: bool,
+
+    /// Writes the matrix files and stop
+    #[structopt(long)]
+    write_matrix: bool,
+
+    /// Ordering strategy
+    #[structopt(short = "o", long, default_value = "Auto")]
+    ordering: String,
+
+    /// Scaling strategy
+    #[structopt(short = "s", long, default_value = "Auto")]
+    scaling: String,
 }
 
 fn main() -> Result<(), StrError> {
@@ -56,10 +68,15 @@ fn main() -> Result<(), StrError> {
     // set the number of BLAS threads
     set_num_threads(opt.blas_nt);
 
-    // get get ODE system
+    // ODE system
     let alpha = if opt.first_book { 2e-3 } else { 0.1 };
     let (system, mut data, mut args) = Samples::brusselator_pde(alpha, opt.npoint, !opt.first_book, false);
 
+    // parameters for the linear solver
+    let mut ls_params = LinSolParams::new();
+    ls_params.ordering = Ordering::from(&opt.ordering);
+    ls_params.scaling = Scaling::from(&opt.scaling);
+
     // set configuration parameters
     let mut params = if opt.dopri8 {
         Params::new(Method::DoPri8)
@@ -73,6 +90,10 @@ fn main() -> Result<(), StrError> {
     params.radau5.concurrent = !opt.serial;
     params.set_tolerances(tol, tol, None)?;
     params.newton.genie = if opt.mumps { Genie::Mumps } else { Genie::Umfpack };
+    params.newton.lin_sol_params = Some(ls_params);
+    if opt.write_matrix {
+        params.newton.write_matrix_after_nstep_and_stop = Some(0);
+    }
 
     // solve the ODE system
     let mut solver = OdeSolver::new(params, &system)?;
diff --git a/russell_ode/src/params.rs b/russell_ode/src/params.rs
index 0f595015..287524a0 100644
--- a/russell_ode/src/params.rs
+++ b/russell_ode/src/params.rs
@@ -32,6 +32,25 @@ pub struct ParamsNewton {
 
     /// Configurations for sparse linear solver
     pub lin_sol_params: Option<LinSolParams>,
+
+    /// Writes the Jacobian and coefficient matrices and stop (with an error message)
+    ///
+    /// The files will be written `if n_accepted > nstep`
+    ///
+    /// Will write the following files:
+    ///
+    /// * `/tmp/russell_ode/jacobian.{mtx,smat}` -- Jacobian matrix
+    /// * `/tmp/russell_ode/kk_real.{mtx,smat}` -- Radau5 coefficient matrix of the real system
+    /// * `/tmp/russell_ode/kk_comp.{mtx,smat}` -- Radau5 coefficient matrix of the complex system
+    ///
+    /// where `mtx` is the extension for the MatrixMarket files and `smat` is the extension
+    /// for the vismatrix files (for visualization).
+    ///
+    /// # References
+    ///
+    /// * MatrixMarket: <https://math.nist.gov/MatrixMarket/formats.html>
+    /// * Vismatrix: <https://github.com/cpmech/vismatrix>
+    pub write_matrix_after_nstep_and_stop: Option<usize>,
 }
 
 /// Holds parameters to control the variable stepsize algorithm
@@ -230,6 +249,7 @@ impl ParamsNewton {
             use_numerical_jacobian: false,
             genie: Genie::Umfpack,
             lin_sol_params: None,
+            write_matrix_after_nstep_and_stop: None,
         }
     }
 
diff --git a/russell_ode/src/radau5.rs b/russell_ode/src/radau5.rs
index 39371799..ac7ebf63 100644
--- a/russell_ode/src/radau5.rs
+++ b/russell_ode/src/radau5.rs
@@ -3,7 +3,7 @@ use crate::{OdeSolverTrait, Params, System, Workspace};
 use num_complex::Complex64;
 use russell_lab::math::SQRT_6;
 use russell_lab::{complex_vec_zip, cpx, format_fortran, vec_copy, ComplexVector, Vector};
-use russell_sparse::numerical_jacobian;
+use russell_sparse::{numerical_jacobian, ComplexCscMatrix, CscMatrix};
 use russell_sparse::{ComplexLinSolver, ComplexSparseMatrix, CooMatrix, Genie, LinSolver, SparseMatrix};
 use std::thread;
 
@@ -231,6 +231,22 @@ where
                 }
             }
         }
+
+        // write the matrices and stop
+        if let Some(nstep) = self.params.newton.write_matrix_after_nstep_and_stop {
+            if work.stats.n_accepted > nstep {
+                let csc_jacobian = CscMatrix::from_coo(jj).unwrap();
+                let csc_kk_real = CscMatrix::from_coo(&kk_real).unwrap();
+                let csc_kk_comp = ComplexCscMatrix::from_coo(&kk_comp).unwrap();
+                csc_jacobian.write_matrix_market("/tmp/russell_ode/jacobian.smat", true)?;
+                csc_jacobian.write_matrix_market("/tmp/russell_ode/jacobian.mtx", false)?;
+                csc_kk_real.write_matrix_market("/tmp/russell_ode/kk_real.smat", true)?;
+                csc_kk_real.write_matrix_market("/tmp/russell_ode/kk_real.mtx", false)?;
+                csc_kk_comp.write_matrix_market("/tmp/russell_ode/kk_comp.smat", true)?;
+                csc_kk_comp.write_matrix_market("/tmp/russell_ode/kk_comp.mtx", false)?;
+                return Err("MATRIX FILES GENERATED in /tmp/russell_ode/");
+            }
+        }
         Ok(())
     }