diff --git a/dynare++/tl/testing/factory.hh b/dynare++/tl/testing/factory.hh
index 284dc6f28..bb306837a 100644
--- a/dynare++/tl/testing/factory.hh
+++ b/dynare++/tl/testing/factory.hh
@@ -47,37 +47,29 @@ public:
   }
 
   template <class _Ttype, class _Ctype>
-  _Ctype *
+  _Ctype
   makeCont(int r, const IntSequence &nvs, int maxdim)
   {
     int symnum = nvs.size();
-    auto *res = new _Ctype(symnum);
+    _Ctype res(symnum);
     for (int dim = 1; dim <= maxdim; dim++)
       if (symnum == 1)
-        {
-          // full symmetry
-          Symmetry sym{dim};
-          res->insert(make<_Ttype>(r, sym, nvs));
-        }
+        // full symmetry
+        res.insert(make<_Ttype>(r, Symmetry{dim}, nvs));
       else
-        {
-          // general symmetry
-          for (int i = 0; i <= dim; i++)
-            {
-              Symmetry sym{i, dim-i};
-              res->insert(make<_Ttype>(r, sym, nvs));
-            }
-        }
+        // general symmetry
+        for (int i = 0; i <= dim; i++)
+          res.insert(make<_Ttype>(r, Symmetry{i, dim-i}, nvs));
     return res;
   }
 
   template <class _Ttype, class _Ptype>
-  _Ptype *
+  _Ptype
   makePoly(int r, int nv, int maxdim)
   {
-    auto *p = new _Ptype(r, nv);
+    _Ptype p(r, nv);
     for (int d = 1; d <= maxdim; d++)
-      p->insert(make<_Ttype>(r, nv, d));
+      p.insert(make<_Ttype>(r, nv, d));
     return p;
   }
 
diff --git a/dynare++/tl/testing/monoms.cc b/dynare++/tl/testing/monoms.cc
index 89c33308f..f3951c5a2 100644
--- a/dynare++/tl/testing/monoms.cc
+++ b/dynare++/tl/testing/monoms.cc
@@ -106,7 +106,7 @@ std::unique_ptr<FGSTensor>
 Monom1Vector::deriv(int dim) const
 {
   auto res = std::make_unique<FGSTensor>(len, TensorDimens(Symmetry{dim},
-                                                           IntSequence(1, nx)));
+                                                           IntSequence{nx}));
   for (Tensor::index it = res->begin(); it != res->end(); ++it)
     {
       Vector outcol{res->getCol(*it)};
@@ -189,16 +189,15 @@ Monom2Vector::deriv(const Symmetry &s) const
   return t;
 }
 
-FGSContainer *
+FGSContainer
 Monom2Vector::deriv(int maxdim) const
 {
-  auto *res = new FGSContainer(2);
+  FGSContainer res(2);
   for (int dim = 1; dim <= maxdim; dim++)
     for (int ydim = 0; ydim <= dim; ydim++)
       {
         int udim = dim - ydim;
-        Symmetry s{ydim, udim};
-        res->insert(deriv(s));
+        res.insert(deriv(Symmetry{ydim, udim}));
       }
   return res;
 }
@@ -331,12 +330,12 @@ Monom4Vector::deriv(const Symmetry &s) const
   return res;
 }
 
-FSSparseTensor *
+FSSparseTensor
 Monom4Vector::deriv(int dim) const
 {
   IntSequence cum{0, nx1, nx1+nx2, nx1+nx2+nx3};
 
-  auto *res = new FSSparseTensor(dim, nx1+nx2+nx3+nx4, len);
+  FSSparseTensor res(dim, nx1+nx2+nx3+nx4, len);
 
   FFSTensor dummy(0, nx1+nx2+nx3+nx4, dim);
   for (Tensor::index run = dummy.begin(); run != dummy.end(); ++run)
@@ -356,7 +355,7 @@ Monom4Vector::deriv(int dim) const
       deriv(ind_sym, ind, col);
       for (int i = 0; i < len; i++)
         if (col[i] != 0.0)
-          res->insert(run.getCoor(), i, col[i]);
+          res.insert(run.getCoor(), i, col[i]);
     }
 
   return res;
@@ -384,7 +383,7 @@ Monom4Vector::print() const
 
 SparseDerivGenerator::SparseDerivGenerator(int nf, int ny, int nu, int nup, int nbigg, int ng,
                                            int mx, double prob, int maxdim)
-  : maxdimen(maxdim), bigg(4), g(4), rcont(4), ts(new FSSparseTensor *[maxdimen])
+  : maxdimen(maxdim), bigg(4), g(4), rcont(4)
 {
   intgen.init(nf, ny, nu, nup, nbigg, mx, prob);
 
@@ -403,20 +402,13 @@ SparseDerivGenerator::SparseDerivGenerator(int nf, int ny, int nu, int nup, int
             g.insert(g_m.deriv(si));
         }
 
-      ts[dim-1] = f.deriv(dim);
+      ts.push_back(f.deriv(dim));
     }
 }
 
-SparseDerivGenerator::~SparseDerivGenerator()
-{
-  for (int i = 0; i < maxdimen; i++)
-    delete ts[i];
-  delete [] ts;
-}
-
 DenseDerivGenerator::DenseDerivGenerator(int ng, int nx, int ny, int nu,
                                          int mx, double prob, int maxdim)
-  : maxdimen(maxdim), ts(maxdimen), uts(maxdimen)
+  : maxdimen(maxdim), xcont(0), rcont(0), ts(maxdimen), uxcont(0), uts(maxdimen)
 {
   intgen.init(ng, nx, ny, nu, nu, mx, prob);
   Monom1Vector g(nx, ng);
@@ -424,7 +416,6 @@ DenseDerivGenerator::DenseDerivGenerator(int ng, int nx, int ny, int nu,
   Monom2Vector r(g, x);
   xcont = x.deriv(maxdimen);
   rcont = r.deriv(maxdimen);
-  uxcont = nullptr;
   for (int d = 1; d <= maxdimen; d++)
     ts[d-1] = g.deriv(d);
 }
@@ -432,13 +423,7 @@ DenseDerivGenerator::DenseDerivGenerator(int ng, int nx, int ny, int nu,
 void
 DenseDerivGenerator::unfold()
 {
-  uxcont = new UGSContainer(*xcont);
+  uxcont = UGSContainer(xcont);
   for (int i = 0; i < maxdimen; i++)
     uts[i] = std::make_unique<UGSTensor>(*(ts[i]));
 }
-
-DenseDerivGenerator::~DenseDerivGenerator()
-{
-  delete xcont;
-  delete rcont;
-}
diff --git a/dynare++/tl/testing/monoms.hh b/dynare++/tl/testing/monoms.hh
index b38e4c925..0de3dc7f8 100644
--- a/dynare++/tl/testing/monoms.hh
+++ b/dynare++/tl/testing/monoms.hh
@@ -68,7 +68,7 @@ public:
   ~Monom2Vector() = default;
   void deriv(const Symmetry &s, const IntSequence &c, Vector &out) const;
   std::unique_ptr<FGSTensor> deriv(const Symmetry &s) const;
-  FGSContainer *deriv(int maxdim) const;
+  FGSContainer deriv(int maxdim) const;
   void print() const;
 };
 
@@ -88,7 +88,7 @@ public:
   Monom4Vector(const Monom4Vector &f, const Monom4Vector &bigg,
                const Monom4Vector &g);
   ~Monom4Vector() = default;
-  FSSparseTensor *deriv(int dim) const;
+  FSSparseTensor deriv(int dim) const;
   std::unique_ptr<FGSTensor> deriv(const Symmetry &s) const;
   void deriv(const Symmetry &s, const IntSequence &coor, Vector &out) const;
   void print() const;
@@ -102,24 +102,22 @@ struct SparseDerivGenerator
   FGSContainer bigg;
   FGSContainer g;
   FGSContainer rcont;
-  FSSparseTensor **const ts;
+  std::vector<FSSparseTensor> ts;
   SparseDerivGenerator(int nf, int ny, int nu, int nup, int nbigg, int ng,
                        int mx, double prob, int maxdim);
-  ~SparseDerivGenerator();
 };
 
 struct DenseDerivGenerator
 {
   int maxdimen;
-  FGSContainer *xcont;
-  FGSContainer *rcont;
+  FGSContainer xcont;
+  FGSContainer rcont;
   std::vector<std::unique_ptr<FGSTensor>> ts;
-  UGSContainer *uxcont;
+  UGSContainer uxcont;
   std::vector<std::unique_ptr<UGSTensor>> uts;
   DenseDerivGenerator(int ng, int nx, int ny, int nu,
                       int mx, double prob, int maxdim);
   void unfold();
-  ~DenseDerivGenerator();
 };
 
 #endif
diff --git a/dynare++/tl/testing/tests.cc b/dynare++/tl/testing/tests.cc
index 6d9eaa7d8..d1dbac829 100644
--- a/dynare++/tl/testing/tests.cc
+++ b/dynare++/tl/testing/tests.cc
@@ -113,18 +113,18 @@ TestRunnable::index_forward(const Symmetry &s, const IntSequence &nvs)
   int ndecr = 0;
   int nincr = 0;
   _Ttype dummy(0, TensorDimens(s, nvs));
-  typename _Ttype::index run = dummy.end();
+  auto run = dummy.end();
   do
     {
       --run;
       ndecr++;
-      typename _Ttype::index run2 = dummy.begin();
+      auto run2 = dummy.begin();
       for (int i = 0; i < *run; i++)
         {
           ++run2;
           nincr++;
         }
-      if (!(run == run2))
+      if (run != run2)
         fails++;
     }
   while (run != dummy.begin());
@@ -145,16 +145,16 @@ TestRunnable::index_backward(const Symmetry &s, const IntSequence &nvs)
   int ndecr = 0;
   int nincr = 0;
   _Ttype dummy(0, TensorDimens(s, nvs));
-  typename _Ttype::index run = dummy.begin();
+  auto run = dummy.begin();
   while (run != dummy.end())
     {
-      typename _Ttype::index run2 = dummy.end();
+      auto run2 = dummy.end();
       for (int i = 0; i < dummy.ncols() - *run; i++)
         {
           --run2;
           ndecr++;
         }
-      if (!(run == run2))
+      if (run != run2)
         fails++;
       ++run;
       nincr++;
@@ -175,11 +175,10 @@ TestRunnable::index_offset(const Symmetry &s, const IntSequence &nvs)
   int fails = 0;
   int nincr = 0;
   _Ttype dummy(0, TensorDimens(s, nvs));
-  for (typename _Ttype::index run = dummy.begin();
-       run != dummy.end(); ++run, nincr++)
+  for (auto run = dummy.begin(); run != dummy.end(); ++run, nincr++)
     {
       typename _Ttype::index run2(dummy, run.getCoor());
-      if (!(run == run2))
+      if (run != run2)
         fails++;
     }
 
@@ -211,16 +210,15 @@ TestRunnable::dense_prod(const Symmetry &bsym, const IntSequence &bnvs,
                          int hdim, int hnv, int rows)
 {
   Factory f;
-  FGSContainer *cont
-    = f.makeCont<FGSTensor, FGSContainer>(hnv, bnvs, bsym.dimen()-hdim+1);
+  FGSContainer cont{f.makeCont<FGSTensor, FGSContainer>(hnv, bnvs, bsym.dimen()-hdim+1)};
   auto fh = f.make<FGSTensor>(rows, Symmetry{hdim}, IntSequence(1, hnv));
   UGSTensor uh(*fh);
   FGSTensor fb(rows, TensorDimens(bsym, bnvs));
   fb.getData().zeros();
   clock_t s1 = clock();
-  cont->multAndAdd(uh, fb);
+  cont.multAndAdd(uh, fb);
   clock_t s2 = clock();
-  UGSContainer ucont(*cont);
+  UGSContainer ucont(cont);
   clock_t s3 = clock();
   UGSTensor ub(rows, fb.getDims());
   ub.getData().zeros();
@@ -241,8 +239,6 @@ TestRunnable::dense_prod(const Symmetry &bsym, const IntSequence &bnvs,
             << "\tunfolded difference norm1:   " << norm1 << '\n'
             << "\tunfolded difference normInf: " << normInf << '\n';
 
-  delete cont;
-
   return norm < 1.e-13;
 }
 
@@ -265,11 +261,11 @@ TestRunnable::folded_monomial(int ng, int nx, int ny, int nu, int dim)
       res.getData().zeros();
       clock_t stime = clock();
       for (int d = 1; d <= dim; d++)
-        gen.xcont->multAndAdd(*(gen.ts[d-1]), res);
+        gen.xcont.multAndAdd(*(gen.ts[d-1]), res);
       stime = clock() - stime;
       std::cout << "\t\ttime for symmetry: "
                 << static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
-      const FGSTensor &mres = gen.rcont->get(s);
+      const FGSTensor &mres = gen.rcont.get(s);
       res.add(-1.0, mres);
       double normtmp = res.getData().getMax();
       std::cout << "\t\terror normMax:     " << normtmp << '\n';
@@ -303,11 +299,11 @@ TestRunnable::unfolded_monomial(int ng, int nx, int ny, int nu, int dim)
       res.getData().zeros();
       clock_t stime = clock();
       for (int d = 1; d <= dim; d++)
-        gen.uxcont->multAndAdd(*(gen.uts[d-1]), res);
+        gen.uxcont.multAndAdd(*(gen.uts[d-1]), res);
       stime = clock() - stime;
       std::cout << "\t\ttime for symmetry: "
                 << static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
-      const FGSTensor &mres = gen.rcont->get(s);
+      const FGSTensor &mres = gen.rcont.get(s);
       FGSTensor foldres(res);
       foldres.add(-1.0, mres);
       double normtmp = foldres.getData().getMax();
@@ -329,7 +325,7 @@ TestRunnable::fold_zcont(int nf, int ny, int nu, int nup, int nbigg,
   for (int d = 1; d <= dim; d++)
     std::cout << "\tfill of dim=" << d << " tensor:     "
               << std::setprecision(2) << std::fixed << std::setw(6)
-              << 100*dg.ts[d-1]->getFillFactor()
+              << 100*dg.ts[d-1].getFillFactor()
               << std::setprecision(6) << std::defaultfloat << " %\n";
   std::cout << "\ttime for monom generation: "
             << static_cast<double>(gen_time)/CLOCKS_PER_SEC << '\n';
@@ -349,7 +345,7 @@ TestRunnable::fold_zcont(int nf, int ny, int nu, int nup, int nbigg,
         res.getData().zeros();
         clock_t stime = clock();
         for (int l = 1; l <= si.dimen(); l++)
-          zc.multAndAdd(*(dg.ts[l-1]), res);
+          zc.multAndAdd(dg.ts[l-1], res);
         stime = clock() - stime;
         std::cout << "\t\ttime for symmetry: "
                   << static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
@@ -375,7 +371,7 @@ TestRunnable::unfold_zcont(int nf, int ny, int nu, int nup, int nbigg,
   for (int d = 1; d <= dim; d++)
     std::cout << "\tfill of dim=" << d << " tensor:     "
               << std::setprecision(2) << std::fixed << std::setw(6)
-              << 100*dg.ts[d-1]->getFillFactor()
+              << 100*dg.ts[d-1].getFillFactor()
               << std::setprecision(6) << std::defaultfloat << " %\n";
   std::cout << "\ttime for monom generation: "
             << static_cast<double>(gen_time)/CLOCKS_PER_SEC << '\n';
@@ -402,7 +398,7 @@ TestRunnable::unfold_zcont(int nf, int ny, int nu, int nup, int nbigg,
         res.getData().zeros();
         clock_t stime = clock();
         for (int l = 1; l <= si.dimen(); l++)
-          zc.multAndAdd(*(dg.ts[l-1]), res);
+          zc.multAndAdd(dg.ts[l-1], res);
         stime = clock() - stime;
         std::cout << "\t\ttime for symmetry: "
                   << static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
@@ -425,14 +421,10 @@ TestRunnable::folded_contraction(int r, int nv, int dim)
   Vector x{fact.makeVector(nv)};
 
   auto forig = fact.make<FFSTensor>(r, nv, dim);
-  auto *f = new FFSTensor(*forig);
+  auto f = std::make_unique<FFSTensor>(*forig);
   clock_t ctime = clock();
   for (int d = dim-1; d > 0; d--)
-    {
-      FFSTensor *fnew = new FFSTensor(*f, ConstVector(x));
-      delete f;
-      f = fnew;
-    }
+    f = std::make_unique<FFSTensor>(*f, ConstVector(x));
   ctime = clock() - ctime;
   Vector res(forig->nrows());
   res.zeros();
@@ -454,8 +446,6 @@ TestRunnable::folded_contraction(int r, int nv, int dim)
             << "\terror normMax:     " << v.getMax() << '\n'
             << "\terror norm1:       " << v.getNorm1() << '\n';
 
-  delete f;
-
   return (v.getMax() < 1.e-10);
 }
 
@@ -467,14 +457,10 @@ TestRunnable::unfolded_contraction(int r, int nv, int dim)
 
   auto forig = fact.make<FFSTensor>(r, nv, dim);
   UFSTensor uorig(*forig);
-  auto *u = new UFSTensor(uorig);
+  auto u = std::make_unique<UFSTensor>(uorig);
   clock_t ctime = clock();
   for (int d = dim-1; d > 0; d--)
-    {
-      UFSTensor *unew = new UFSTensor(*u, ConstVector(x));
-      delete u;
-      u = unew;
-    }
+    u = std::make_unique<UFSTensor>(*u, ConstVector(x));
   ctime = clock() - ctime;
   Vector res(uorig.nrows());
   res.zeros();
@@ -495,8 +481,6 @@ TestRunnable::unfolded_contraction(int r, int nv, int dim)
             << "\terror normMax:     " << v.getMax() << '\n'
             << "\terror norm1:       " << v.getNorm1() << '\n';
 
-  delete u;
-
   return (v.getMax() < 1.e-10);
 }
 
@@ -515,34 +499,30 @@ TestRunnable::poly_eval(int r, int nv, int maxdim)
   Vector out_uh(r);
   out_uh.zeros();
 
-  UTensorPolynomial *up;
-  {
-    FTensorPolynomial *fp = fact.makePoly<FFSTensor, FTensorPolynomial>(r, nv, maxdim);
+  FTensorPolynomial fp{fact.makePoly<FFSTensor, FTensorPolynomial>(r, nv, maxdim)};
 
-    clock_t ft_cl = clock();
-    fp->evalTrad(out_ft, x);
-    ft_cl = clock() - ft_cl;
-    std::cout << "\ttime for folded power eval:    "
-              << static_cast<double>(ft_cl)/CLOCKS_PER_SEC << '\n';
+  clock_t ft_cl = clock();
+  fp.evalTrad(out_ft, x);
+  ft_cl = clock() - ft_cl;
+  std::cout << "\ttime for folded power eval:    "
+    << static_cast<double>(ft_cl)/CLOCKS_PER_SEC << '\n';
 
-    clock_t fh_cl = clock();
-    fp->evalHorner(out_fh, x);
-    fh_cl = clock() - fh_cl;
-    std::cout << "\ttime for folded horner eval:   "
-              << static_cast<double>(fh_cl)/CLOCKS_PER_SEC << '\n';
+  clock_t fh_cl = clock();
+  fp.evalHorner(out_fh, x);
+  fh_cl = clock() - fh_cl;
+  std::cout << "\ttime for folded horner eval:   "
+    << static_cast<double>(fh_cl)/CLOCKS_PER_SEC << '\n';
 
-    up = new UTensorPolynomial(*fp);
-    delete fp;
-  }
+  UTensorPolynomial up{fp};
 
   clock_t ut_cl = clock();
-  up->evalTrad(out_ut, x);
+  up.evalTrad(out_ut, x);
   ut_cl = clock() - ut_cl;
   std::cout << "\ttime for unfolded power eval:  "
             << static_cast<double>(ut_cl)/CLOCKS_PER_SEC << '\n';
 
   clock_t uh_cl = clock();
-  up->evalHorner(out_uh, x);
+  up.evalHorner(out_uh, x);
   uh_cl = clock() - uh_cl;
   std::cout << "\ttime for unfolded horner eval: "
             << static_cast<double>(uh_cl)/CLOCKS_PER_SEC << '\n';
@@ -558,7 +538,6 @@ TestRunnable::poly_eval(int r, int nv, int maxdim)
             << "\tfolded horner error norm max:    " << max_fh << '\n'
             << "\tunfolded horner error norm max:  " << max_uh << '\n';
 
-  delete up;
   return (max_ft+max_fh+max_uh < 1.0e-10);
 }