diff --git a/SpMat.h b/SpMat.h
index 0b35796b25d801fbf15f80d310a1d07996b17fa9..07377c7b19f38ff9fb82750a8dc3d8db653f405a 100644
--- a/SpMat.h
+++ b/SpMat.h
@@ -21,10 +21,10 @@
 #ifndef SpMat_h
 #define SpMat_h
 
+#include <memory>
 #include <vector>
 #include <fstream>
 #include <iomanip>
-#include <boost/shared_ptr.hpp>
 #include "armawrap/newmat.h"
 #include "cg.h"
 #include "bicg.h"
@@ -188,7 +188,7 @@ public:
                                  MatrixType                            type = UNKNOWN,
                                  double                                tol = 1e-4,
                                  unsigned int                          miter = 200,
-				 boost::shared_ptr<Preconditioner<T> > C = boost::shared_ptr<Preconditioner<T> >()) const;
+				 std::shared_ptr<Preconditioner<T> > C = std::shared_ptr<Preconditioner<T> >()) const;
 
   NEWMAT::ReturnMatrix SolveForx(const NEWMAT::ColumnVector&            b,
 			         MatrixType                             type,
@@ -200,7 +200,7 @@ public:
 			         MatrixType                             type,
 			         double                                 tol,
 			         unsigned int                           miter,
-				 boost::shared_ptr<Preconditioner<T> >  C,
+				 std::shared_ptr<Preconditioner<T> >  C,
 				 const NEWMAT::ColumnVector&            x_init) const;
 
 protected:
@@ -559,7 +559,7 @@ NEWMAT::ReturnMatrix SpMat<T>::SolveForx(const NEWMAT::ColumnVector&
 			                 unsigned int                           miter,
 					 const NEWMAT::ColumnVector&            x_init) const
 {
-  return this->SolveForx(b,type,tol,miter,boost::shared_ptr<Preconditioner<T> >(),x_init);
+  return this->SolveForx(b,type,tol,miter,std::shared_ptr<Preconditioner<T> >(),x_init);
 }
 
 template<class T>
@@ -567,7 +567,7 @@ NEWMAT::ReturnMatrix SpMat<T>::SolveForx(const NEWMAT::ColumnVector&
 			                 MatrixType                             type,
 			                 double                                 tol,
 			                 unsigned int                           miter,
-					 boost::shared_ptr<Preconditioner<T> >  C) const
+					 std::shared_ptr<Preconditioner<T> >  C) const
 {
   NEWMAT::ColumnVector x_init;
   return this->SolveForx(b,type,tol,miter,C,x_init);
@@ -578,7 +578,7 @@ NEWMAT::ReturnMatrix SpMat<T>::SolveForx(const NEWMAT::ColumnVector&
 			                 MatrixType                             type,
 			                 double                                 tol,
 			                 unsigned int                           miter,
-					 boost::shared_ptr<Preconditioner<T> >  C,
+					 std::shared_ptr<Preconditioner<T> >  C,
 					 const NEWMAT::ColumnVector&           x_init) const
 {
   if (_m != _n) throw SpMatException("SolveForx: Matrix must be square");
@@ -598,8 +598,8 @@ NEWMAT::ReturnMatrix SpMat<T>::SolveForx(const NEWMAT::ColumnVector&
   int                   liter = int(miter);
   double                ltol = tol;
   // Use diagonal conditioner if no user-specified one
-  boost::shared_ptr<Preconditioner<T> > M = boost::shared_ptr<Preconditioner<T> >();
-  if (!C) M = boost::shared_ptr<Preconditioner<T> >(new DiagPrecond<T>(*this));
+  std::shared_ptr<Preconditioner<T> > M = std::shared_ptr<Preconditioner<T> >();
+  if (!C) M = std::shared_ptr<Preconditioner<T> >(new DiagPrecond<T>(*this));
   else M = C;
 
   switch (type) {
diff --git a/bfmatrix.cpp b/bfmatrix.cpp
index 81bd9689a0c0c8096964a4655cbea1cafa492516..d8f34efdebf7500504a62ae347dcc19504de2685 100644
--- a/bfmatrix.cpp
+++ b/bfmatrix.cpp
@@ -6,12 +6,11 @@
 //    Copyright (C) 2007 University of Oxford
 //
 /*  CCOPYRIGHT  */
+#include <memory>
 #include <iostream>
 #include <iomanip>
 #include <fstream>
-#include <boost/shared_ptr.hpp>
 #include "armawrap/newmat.h"
-#include "armawrap/newmatio.h"
 #include "miscmaths.h"
 #include "bfmatrix.h"
 
@@ -45,10 +44,10 @@ void FullBFMatrix::Print(const std::string fname) const
   else write_ascii_matrix(fname,*mp);
 }
 
-boost::shared_ptr<BFMatrix> FullBFMatrix::Transpose()
+std::shared_ptr<BFMatrix> FullBFMatrix::Transpose()
 const
 {
-  boost::shared_ptr<FullBFMatrix>  tm(new FullBFMatrix(mp->t()));
+  std::shared_ptr<FullBFMatrix>  tm(new FullBFMatrix(mp->t()));
   return(tm);
 }
 
diff --git a/bfmatrix.h b/bfmatrix.h
index 673bfc82653e60948e2b324a6f8e64eaf5966248..2ba2d7f6ec544a1fa948f43cf2a0b566097fe531 100644
--- a/bfmatrix.h
+++ b/bfmatrix.h
@@ -28,7 +28,7 @@
 #ifndef BFMatrix_h
 #define BFMatrix_h
 
-#include <boost/shared_ptr.hpp>
+#include <memory>
 #include "armawrap/newmat.h"
 #include "SpMat.h"
 #include "cg.h"
@@ -98,7 +98,7 @@ public:
   virtual void AddTo(unsigned int x, unsigned int y, double val) = 0;
 
   // Transpose
-  virtual boost::shared_ptr<BFMatrix> Transpose() const = 0;
+  virtual std::shared_ptr<BFMatrix> Transpose() const = 0;
 
   // Concatenation. Note that desired polymorphism prevents us from using BFMatrix->NEWMAT::Matrix conversion
   // Concatenate two matrices yielding a third
@@ -131,25 +131,25 @@ template<class T>
 class SparseBFMatrix : public BFMatrix
 {
 private:
-  boost::shared_ptr<MISCMATHS::SpMat<T> >    mp;
+  std::shared_ptr<MISCMATHS::SpMat<T> >    mp;
 
 public:
   // Constructors, destructor and assignment
   SparseBFMatrix()
-  : mp(boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>())) {}
+  : mp(std::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>())) {}
   SparseBFMatrix(unsigned int m, unsigned int n)
-  : mp(boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(m,n))) {}
+  : mp(std::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(m,n))) {}
   SparseBFMatrix(unsigned int m, unsigned int n, const unsigned int *irp, const unsigned int *jcp, const double *sp)
-  : mp(boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(m,n,irp,jcp,sp))) {}
+  : mp(std::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(m,n,irp,jcp,sp))) {}
   SparseBFMatrix(const MISCMATHS::SpMat<T>& M)
-  : mp(boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(M))) {}
+  : mp(std::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(M))) {}
   SparseBFMatrix(const NEWMAT::Matrix& M)
-  : mp(boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(M))) {}
+  : mp(std::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(M))) {}
   SparseBFMatrix(const SparseBFMatrix<T>& M)
-  : mp(boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(*(M.mp)))) {}
+  : mp(std::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(*(M.mp)))) {}
   virtual ~SparseBFMatrix() {}
   virtual const SparseBFMatrix& operator=(const SparseBFMatrix<T>& M) {
-    mp = boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(*(M.mp))); return(*this);
+    mp = std::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(*(M.mp))); return(*this);
   }
 
   friend class BFMatrixColumnIterator;
@@ -165,12 +165,12 @@ public:
   virtual void Print(const std::string fname=std::string("")) const {mp->Print(fname);}
 
   // Setting, deleting or resizing the whole sparse matrix.
-  virtual void SetMatrix(const MISCMATHS::SpMat<T>& M) {mp = boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(M));}
-  // virtual void SetMatrix(const MISCMATHS::SpMat<float>& M) {mp = boost::shared_ptr<MISCMATHS::SpMat<float> >(new MISCMATHS::SpMat<float>(M));}
-  virtual void SetMatrix(const NEWMAT::Matrix& M) {mp = boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(M));}
-  virtual void SetMatrixPtr(boost::shared_ptr<MISCMATHS::SpMat<T> >& mptr) {mp = mptr;}
-  virtual void Clear() {mp = boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>());}
-  virtual void Resize(unsigned int m, unsigned int n) {mp = boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(m,n));}
+  virtual void SetMatrix(const MISCMATHS::SpMat<T>& M) {mp = std::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(M));}
+  // virtual void SetMatrix(const MISCMATHS::SpMat<float>& M) {mp = std::shared_ptr<MISCMATHS::SpMat<float> >(new MISCMATHS::SpMat<float>(M));}
+  virtual void SetMatrix(const NEWMAT::Matrix& M) {mp = std::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(M));}
+  virtual void SetMatrixPtr(std::shared_ptr<MISCMATHS::SpMat<T> >& mptr) {mp = mptr;}
+  virtual void Clear() {mp = std::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>());}
+  virtual void Resize(unsigned int m, unsigned int n) {mp = std::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(m,n));}
 
   // Accessing values
   virtual double Peek(unsigned int r, unsigned int c) const {return(mp->Peek(r,c));}
@@ -181,7 +181,7 @@ public:
   virtual void AddTo(unsigned int x, unsigned int y, double val) {mp->AddTo(x,y,val);}
 
   // Transpose.
-  virtual boost::shared_ptr<BFMatrix> Transpose() const;
+  virtual std::shared_ptr<BFMatrix> Transpose() const;
 
   // Concatenation of two matrices returning a third
   // AB = [*this B] in Matlab lingo
@@ -216,16 +216,16 @@ public:
 class FullBFMatrix : public BFMatrix
 {
 private:
-  boost::shared_ptr<NEWMAT::Matrix>    mp;
+  std::shared_ptr<NEWMAT::Matrix>    mp;
 public:
   // Constructors, destructor and assignment
-  FullBFMatrix() {mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix());}
-  FullBFMatrix(unsigned int m, unsigned int n) {mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(m,n));}
-  FullBFMatrix(const MISCMATHS::SpMat<double>& M) {mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M.AsNEWMAT()));}
-  FullBFMatrix(const NEWMAT::Matrix& M) {mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M));}
+  FullBFMatrix() {mp = std::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix());}
+  FullBFMatrix(unsigned int m, unsigned int n) {mp = std::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(m,n));}
+  FullBFMatrix(const MISCMATHS::SpMat<double>& M) {mp = std::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M.AsNEWMAT()));}
+  FullBFMatrix(const NEWMAT::Matrix& M) {mp = std::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M));}
   virtual ~FullBFMatrix() {}
   virtual const FullBFMatrix& operator=(const FullBFMatrix& M) {
-    mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(*(M.mp))); return(*this);
+    mp = std::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(*(M.mp))); return(*this);
   }
 
   friend class BFMatrixColumnIterator;
@@ -241,10 +241,10 @@ public:
   virtual void Print(const std::string fname=std::string("")) const;
 
   // Setting, deleting or resizing the whole matrix.
-  virtual void SetMatrix(const MISCMATHS::SpMat<double>& M) {mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M.AsNEWMAT()));}
-  virtual void SetMatrix(const MISCMATHS::SpMat<float>& M) {mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M.AsNEWMAT()));}
-  virtual void SetMatrix(const NEWMAT::Matrix& M) {mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M));}
-  virtual void SetMatrixPtr(boost::shared_ptr<NEWMAT::Matrix>& mptr) {mp = mptr;}
+  virtual void SetMatrix(const MISCMATHS::SpMat<double>& M) {mp = std::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M.AsNEWMAT()));}
+  virtual void SetMatrix(const MISCMATHS::SpMat<float>& M) {mp = std::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M.AsNEWMAT()));}
+  virtual void SetMatrix(const NEWMAT::Matrix& M) {mp = std::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M));}
+  virtual void SetMatrixPtr(std::shared_ptr<NEWMAT::Matrix>& mptr) {mp = mptr;}
   virtual void Clear() {mp->ReSize(0,0);}
   virtual void Resize(unsigned int m, unsigned int n) {mp->ReSize(m,n);}
 
@@ -257,7 +257,7 @@ public:
   virtual void AddTo(unsigned int x, unsigned int y, double val) {(*mp)(x,y)+=val;}
 
   // Transpose.
-  virtual boost::shared_ptr<BFMatrix> Transpose() const;
+  virtual std::shared_ptr<BFMatrix> Transpose() const;
 
   // Concatenation of two matrices returning a third
   virtual void HorConcat(const BFMatrix& B, BFMatrix& AB) const;
@@ -383,9 +383,9 @@ private:
 //
 
 template<class T>
-boost::shared_ptr<BFMatrix> SparseBFMatrix<T>::Transpose() const
+std::shared_ptr<BFMatrix> SparseBFMatrix<T>::Transpose() const
 {
-  boost::shared_ptr<SparseBFMatrix<T> >   tm(new SparseBFMatrix<T>(mp->t()));
+  std::shared_ptr<SparseBFMatrix<T> >   tm(new SparseBFMatrix<T>(mp->t()));
   return(tm);
 }
 
diff --git a/nonlin.cpp b/nonlin.cpp
index 8e14cb4df8f591722a19b3993f5952863dbf3872..24fefd0313b943e3e0434b0797647a6fb4089a04 100644
--- a/nonlin.cpp
+++ b/nonlin.cpp
@@ -3,6 +3,7 @@
 /*  CCOPYRIGHT  */
 // Definitions for module nonlin
 
+#include <memory>
 #include <ctime>
 #include <iostream>
 #include <fstream>
@@ -12,10 +13,10 @@
 #include <cmath>
 #include "armawrap/newmat.h"
 #include "armawrap/newmatio.h"
+#include "utils/fsl_isfinite.h"
 #include "bfmatrix.h"
 #include "nonlin.h"
 #include "Simplex.h"
-#include "utils/fsl_isfinite.h"
 
 using namespace std;
 using namespace NEWMAT;
@@ -168,12 +169,12 @@ ReturnMatrix NonlinCF::grad(const ColumnVector& p) const
 // used interchangeably depending on if the structure of the Hessian
 // renders it sparse or not.
 
-boost::shared_ptr<BFMatrix> NonlinCF::hess(const ColumnVector&         p,
-                                           boost::shared_ptr<BFMatrix> iptr) const
+std::shared_ptr<BFMatrix> NonlinCF::hess(const ColumnVector&         p,
+                                           std::shared_ptr<BFMatrix> iptr) const
 {
-  boost::shared_ptr<BFMatrix>    hessm;
+  std::shared_ptr<BFMatrix>    hessm;
   if (iptr && int(iptr->Nrows())==p.Nrows() && int(iptr->Ncols())==p.Nrows()) hessm = iptr;
-  else hessm = boost::shared_ptr<BFMatrix>(new FullBFMatrix(p.Nrows(),p.Nrows()));
+  else hessm = std::shared_ptr<BFMatrix>(new FullBFMatrix(p.Nrows(),p.Nrows()));
   ColumnVector tmpp = p;
   double tiny = 1e-4;
   double fx0y0 = cf(tmpp);
@@ -343,7 +344,7 @@ NonlinOut levmar(const NonlinParam& p, const NonlinCF& cfo)
   bool                            success = true;             // True if last step decreased CF
   double                          olambda = 0.0;              // How much the diagonal of H was nudged last time
   ColumnVector                    g;                          // Gradient
-  boost::shared_ptr<BFMatrix>     H;                          // Hessian
+  std::shared_ptr<BFMatrix>     H;                          // Hessian
 
   while (p.NextIter(success)) {
     if (success) {                                            // If last attempt decreased cost-function
@@ -358,7 +359,7 @@ NonlinOut levmar(const NonlinParam& p, const NonlinCF& cfo)
       else if (p.GaussNewtonType() == LM_L) {                // If Levenberg
         H->AddTo(i,i,p.Lambda()-olambda);
       }
-      else if (p.GaussNewtonType() == LM_GN) {}              // If it is pure Gauss-Newton we don't nudge at all 
+      else if (p.GaussNewtonType() == LM_GN) {}              // If it is pure Gauss-Newton we don't nudge at all
     }
     ColumnVector step;
     double ncf = 0.0;
@@ -1119,10 +1120,10 @@ pair<ColumnVector,ColumnVector> check_grad(const ColumnVector&  par,
   return(rv);
 }
 
-pair<boost::shared_ptr<BFMatrix>,boost::shared_ptr<BFMatrix> > check_hess(const ColumnVector& par,
+pair<std::shared_ptr<BFMatrix>,std::shared_ptr<BFMatrix> > check_hess(const ColumnVector& par,
                                                                           const NonlinCF&     cfo)
 {
-  pair<boost::shared_ptr<BFMatrix>,boost::shared_ptr<BFMatrix> > rv;
+  pair<std::shared_ptr<BFMatrix>,std::shared_ptr<BFMatrix> > rv;
 
   rv.first = cfo.NonlinCF::hess(par);
   rv.second = cfo.hess(par);
diff --git a/nonlin.h b/nonlin.h
index cd23bbe305d93f07798a03a094f3a26336ef66ee..a0d0fdfa658471cba04212a932dbeec5c1057b84 100644
--- a/nonlin.h
+++ b/nonlin.h
@@ -16,11 +16,11 @@
 #ifndef nonlin_h
 #define nonlin_h
 
+#include <memory>
 #include <string>
 #include <vector>
 #include <iostream>
 #include <iomanip>
-#include <boost/shared_ptr.hpp>
 #include "bfmatrix.h"
 #include "armawrap/newmat.h"
 #include "armawrap/newmatio.h"
@@ -328,8 +328,8 @@ public:
   virtual ~NonlinCF() {}
   virtual double sf() const {return(1.0);}
   virtual NEWMAT::ReturnMatrix grad(const NEWMAT::ColumnVector& p) const;
-  virtual boost::shared_ptr<BFMatrix> hess(const NEWMAT::ColumnVector& p,
-                                           boost::shared_ptr<BFMatrix> iptr=boost::shared_ptr<BFMatrix>()) const;
+  virtual std::shared_ptr<BFMatrix> hess(const NEWMAT::ColumnVector& p,
+                                           std::shared_ptr<BFMatrix> iptr=std::shared_ptr<BFMatrix>()) const;
   virtual double cf(const NEWMAT::ColumnVector& p) const = 0;
 };
 
@@ -400,7 +400,7 @@ NonlinOut nonlin(const NonlinParam& p, const NonlinCF& cfo);
 std::pair<NEWMAT::ColumnVector,NEWMAT::ColumnVector> check_grad(const NEWMAT::ColumnVector&  par,
                                                                 const NonlinCF&              cfo);
 
-std::pair<boost::shared_ptr<BFMatrix>,boost::shared_ptr<BFMatrix> > check_hess(const NEWMAT::ColumnVector& par,
+std::pair<std::shared_ptr<BFMatrix>,std::shared_ptr<BFMatrix> > check_hess(const NEWMAT::ColumnVector& par,
                                                                                const NonlinCF&             cfo);
 
 } // End namespace MISCMATHS