diff --git a/src/algebra/LinearEquation.cpp b/src/algebra/LinearEquation.cpp
index afaface63..fa402fc6c 100644
--- a/src/algebra/LinearEquation.cpp
+++ b/src/algebra/LinearEquation.cpp
@@ -42,8 +42,6 @@ LinearEquation::LinearEquation(Mesh *other_msh){
   _EPSC = NULL;
   _RES = NULL;
   _RESC = NULL;
-  _PP = NULL;
-  _RR = NULL;
   _KK = NULL;
   _CC = NULL;
 }
@@ -237,14 +235,6 @@ void LinearEquation::DeletePde() {
       delete _CC;
   }
   
-  
-  if (_msh->GetLevel()>0) {
-     if(_PP) 
-       delete _PP;
-     if(_RR) 
-       delete _RR;
-  }
-  
   if(_EPS)
     delete _EPS;
   
diff --git a/src/algebra/LinearEquation.hpp b/src/algebra/LinearEquation.hpp
index 50bc24a49..8576faad3 100644
--- a/src/algebra/LinearEquation.hpp
+++ b/src/algebra/LinearEquation.hpp
@@ -79,7 +79,7 @@ class LinearEquation : public ParallelObject {
   // member data
   Mesh *_msh; 
   NumericVector *_EPS, *_EPSC, *_RES, *_RESC;
-  SparseMatrix *_KK, *_PP,*_RR, *_CC; 
+  SparseMatrix *_KK, *_CC; 
   vector < vector <unsigned> > KKoffset;
   vector < unsigned > KKghostsize;
   vector < vector < int> > KKghost_nd;
diff --git a/src/equations/LinearImplicitSystem.cpp b/src/equations/LinearImplicitSystem.cpp
index b3b5e4351..55bf29a6d 100644
--- a/src/equations/LinearImplicitSystem.cpp
+++ b/src/equations/LinearImplicitSystem.cpp
@@ -109,6 +109,8 @@ void LinearImplicitSystem::init() {
     AddVariableToBeSolved("All");
 }
 
+
+
 /// @deprecated 
 // this function is like init but it doesn't call InitPDE
 void LinearImplicitSystem::init_two() {
@@ -124,6 +126,15 @@ void LinearImplicitSystem::init_two() {
 //       _LinSolver[i]->InitPde(_SolSystemPdeIndex,_ml_sol->GetSolType(),
 // 			     _ml_sol->GetSolName(),&_solution[i]->_Bdc,_gridr,_gridn,_SparsityPattern);
 //     }  
+
+
+    _PP.resize(_gridn);
+    _RR.resize(_gridn);
+    for(unsigned i=0;i<_gridn;i++){
+      _PP[i]=NULL;
+      _RR[i]=NULL;
+    }
+    
 //     
 //     for (unsigned ig=1; ig<_gridn; ig++) {
 //       BuildProlongatorMatrix(ig);
diff --git a/src/equations/LinearImplicitSystem.hpp b/src/equations/LinearImplicitSystem.hpp
index af5ea8fd5..9afea01f2 100644
--- a/src/equations/LinearImplicitSystem.hpp
+++ b/src/equations/LinearImplicitSystem.hpp
@@ -168,9 +168,10 @@ class LinearImplicitSystem : public ImplicitSystem {
      /** enforce sparcity pattern for setting uncoupled variables and save on memory allocation **/
     void SetSparsityPattern(vector < bool > other_sparcity_pattern);
     
+    vector < SparseMatrix* > _PP, _RR; /// @todo put it back to protected
+
 protected:
 
-    vector < SparseMatrix* > _PP, _RR; 
     /** Create the Prolongator matrix for the Multigrid solver */
     void Prolongator(const unsigned &gridf);
 
diff --git a/src/meshGencase/GenCase.cpp b/src/meshGencase/GenCase.cpp
index 60d437b17..9de00bc4e 100644
--- a/src/meshGencase/GenCase.cpp
+++ b/src/meshGencase/GenCase.cpp
@@ -2115,7 +2115,7 @@ void GenCase::ComputeMaxElXNode() {
 
 
 // =========================================
-void GenCase::ReadMGOps(const std::string output_path, const SystemTwo * mysys) {
+void GenCase::ReadMGOps(const std::string output_path, SystemTwo * mysys) {
 
     std::string     f_matrix = DEFAULT_F_MATRIX;
     std::string       f_rest = DEFAULT_F_REST;
@@ -2297,7 +2297,7 @@ void GenCase::ReadMGOps(const std::string output_path, const SystemTwo * mysys)
 // So, if the node_dof was not filled correctly, then when you 
 
 
-void GenCase::ReadMatrix(const  std::string& namefile, const SystemTwo * mysys) {
+void GenCase::ReadMatrix(const  std::string& namefile, SystemTwo * mysys) {
 
     for (uint Level = 0; Level< mysys->GetGridn(); Level++) {
 
@@ -2490,7 +2490,7 @@ void GenCase::ReadMatrix(const  std::string& namefile, const SystemTwo * mysys)
 
 //=============================
 //This function depends on _iproc
-void GenCase::ReadProl(const std::string& name, const SystemTwo * mysys) {
+void GenCase::ReadProl(const std::string& name, SystemTwo * mysys) {
 
 
     for (uint Level = 1; Level< mysys->GetGridn(); Level++) {
@@ -2573,7 +2573,7 @@ void GenCase::ReadProl(const std::string& name, const SystemTwo * mysys) {
    
     uint off_proc = mysys->GetMLProb().GetMeshTwo()._iproc*mysys->GetGridn();
 
-    mysys->_LinSolver[Lev_f]->_PP = SparseMatrix::build().release();
+    mysys->_PP[Lev_f] = SparseMatrix::build().release();
 // // //     _Prl[ Lev_f ]->init(0,0,0,0); //TODO BACK TO A REASONABLE INIT
 
     // local matrix dimension
@@ -2650,7 +2650,7 @@ void GenCase::ReadProl(const std::string& name, const SystemTwo * mysys) {
 
     std::cout << "Printing Prolongator ===========" << std::endl;
     pattern.print();
-    mysys->_LinSolver[Lev_f]->_PP->update_sparsity_pattern_old(pattern);
+    mysys->_PP[Lev_f]->update_sparsity_pattern_old(pattern);
 
 //=========== VALUES ===================
     DenseMatrix *valmat;
@@ -2670,7 +2670,7 @@ void GenCase::ReadProl(const std::string& name, const SystemTwo * mysys) {
             for (uint j=0; j<ind.size(); j++) ind[j] = pattern[irow][j];
             valmat = new DenseMatrix(1,ncol);
             for (uint j=0; j<ncol; j++)(*valmat)(0,j) = Prol_val[fe][j+len[fe][i]];
-            mysys->_LinSolver[Lev_f]->_PP->add_matrix(*valmat,tmp,ind);
+            mysys->_PP[Lev_f]->add_matrix(*valmat,tmp,ind);
             delete  valmat;
         }
      }
@@ -2691,7 +2691,7 @@ void GenCase::ReadProl(const std::string& name, const SystemTwo * mysys) {
 
     pattern.clear();
 
-    mysys->_LinSolver[Lev_f]->_PP->close();
+    mysys->_PP[Lev_f]->close();
 //     if (mysys->GetMLProb().GetMeshTwo()._iproc==0) _Prl[  Lev_f ]->print_personal();
 //     _Prl[  Lev_f ]->print_graphic(false); //TODO should pass this true or false as a parameter
    } //end levels
@@ -2727,7 +2727,7 @@ void GenCase::ReadProl(const std::string& name, const SystemTwo * mysys) {
     
     //AAA fai molta attenzione: per esplorare la node_dof devi usare Lev_c e Lev_f,
     //perche' sono legati ai DOF (devi pensare che la questione del mesh e' gia' risolta)
-void GenCase::ReadRest(const std::string& name, const SystemTwo * mysys) {
+void GenCase::ReadRest(const std::string& name, SystemTwo * mysys) {
  
   
   for (uint Level = 0; Level< mysys->GetGridn() - 1; Level++) {
@@ -2809,7 +2809,7 @@ void GenCase::ReadRest(const std::string& name, const SystemTwo * mysys) {
 
     uint off_proc=mysys->GetGridn()*mysys->GetMLProb().GetMeshTwo()._iproc;
 
-    mysys->_LinSolver[Lev_c]->_RR = SparseMatrix::build().release();
+    mysys->_RR[Lev_c] = SparseMatrix::build().release();
 // // //     _Rst[Lev_c]->init(0,0,0,0);   //TODO BACK TO A REASONABLE INIT  //we have to do this before appropriately!!!
 
     int nrowt=0;int nclnt=0;
@@ -2873,7 +2873,7 @@ void GenCase::ReadRest(const std::string& name, const SystemTwo * mysys) {
 
     std::cout << "Printing Restrictor ===========" << std::endl;
     pattern.print();
-    mysys->_LinSolver[Lev_c]->_RR->update_sparsity_pattern_old(pattern);  //TODO see 
+    mysys->_RR[Lev_c]->update_sparsity_pattern_old(pattern);  //TODO see 
 //         _Rst[Lev_c]->close();
 //     if (mysys->GetMLProb().GetMeshTwo()._iproc==0) _Rst[Lev_c]->print_personal(); //there is no print function for rectangular matrices, and print_personal doesnt seem to be working...
 // la print stampa il contenuto, ma io voglio solo stampare lo sparsity pattern!
@@ -2901,7 +2901,7 @@ void GenCase::ReadRest(const std::string& name, const SystemTwo * mysys) {
             for (uint i1=0;i1<ind.size();i1++) { ind[i1] = pattern[irow][i1]; /*std::cout << " " << ind[i1] << " ";*/}
             valmat = new DenseMatrix(1,ncol);  //TODO add a matrix row by row...
             for (uint j=0; j<ncol; j++) (*valmat)(0,j) = mysys->_bcond._bc[irow_top]*Rest_val[fe][ j+len[fe][i] ];
-            mysys->_LinSolver[Lev_c]->_RR->add_matrix(*valmat,tmp,ind);
+            mysys->_RR[Lev_c]->add_matrix(*valmat,tmp,ind);
             delete  valmat;
             }// end dof loop
          } // end var loop
@@ -2921,7 +2921,7 @@ void GenCase::ReadRest(const std::string& name, const SystemTwo * mysys) {
 
     pattern.clear();
 
-    mysys->_LinSolver[Lev_c]->_RR->close();
+    mysys->_RR[Lev_c]->close();
 //     if (mysys->GetMLProb().GetMeshTwo()._iproc==0)  _Rst[Lev_c]->print_personal(std::cout);
 //     _Rst[Lev_c]->print_graphic(false); // TODO should pass this true or false as a parameter
 
diff --git a/src/meshGencase/GenCase.hpp b/src/meshGencase/GenCase.hpp
index d7ba36da4..157bf9f21 100644
--- a/src/meshGencase/GenCase.hpp
+++ b/src/meshGencase/GenCase.hpp
@@ -47,10 +47,10 @@ class GenCase : public MultiLevelMeshTwo {
     void ComputeAndPrintProl(const std::string output_path);
     void ComputeAndPrintRest(const std::string output_path);
     
-    static void ReadMGOps(const std::string output_path, const SystemTwo * mysys);
-    static void ReadMatrix(const std::string& name, const SystemTwo * mysys); 
-    static void ReadProl(const std::string& name, const SystemTwo * mysys);   
-    static void ReadRest(const std::string& name, const SystemTwo * mysys);
+    static void ReadMGOps(const std::string output_path, SystemTwo * mysys);
+    static void ReadMatrix(const std::string& name,  SystemTwo * mysys); 
+    static void ReadProl(const std::string& name, SystemTwo * mysys);   
+    static void ReadRest(const std::string& name, SystemTwo * mysys);
     
     void CreateMeshStructuresLevSubd(const std::string output_path);
     void Delete();