--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -15,6 +15,8 @@
 if(CMAKE_BUILD_TYPE STREQUAL "Debug")
   set(PEDANTIC_COMPILE_FLAGS "${PEDANTIC_COMPILE_FLAGS} -pg -g")
   set(LAX_COMPILE_FLAGS "${LAX_COMPILE_FLAGS} -pg")
+  set(GSL_LIBRARIES "gsl_p")             ## profiling 
+  set(GSL_CBLAS_LIBRARIES "gslcblas_p")  ## profiling
 endif()
 
 include_directories("include/")

--- a/src/include/interp_values.hpp
+++ b/src/include/interp_values.hpp
@@ -3,9 +3,11 @@
 
 #include "linalg.hpp"
 #include <map>
+#include <boost/shared_ptr.hpp>
 
 namespace kwantix{
 using std::map;
+using boost::shared_ptr;
 
 class bdry_values;
 
@@ -46,7 +48,7 @@
   size_t m;
 
   ///Store normals in a matrix
-  kwantix::matrix N;
+  shared_ptr<matrix> N;
 };
 
 

--- a/src/include/interpolator.hpp
+++ b/src/include/interpolator.hpp
@@ -215,7 +215,7 @@
   size_t m; 
 
   ///The matrix to invert.
-  matrix M; 
+  shared_ptr<matrix> M; 
 
   ///Is the interpolator ready for use?
   bool initted;
@@ -254,7 +254,8 @@
    * vector represents zeros. Thus, an empty vector represents
    * evaluation instead of derivatives.
    */
-  mutable map<std::vector<size_t>, matrix> precomp_rbfs;
+  mutable map<std::vector<size_t>, 
+              shared_ptr<matrix> > precomp_rbfs;
     
   /** @name  Evaluation flags
    *Have the RBFs or their derivatives been precomputed?

--- a/src/include/linalg.hpp
+++ b/src/include/linalg.hpp
@@ -73,10 +73,16 @@
    *  @param j - Column number.
    *  @return A reference to the matrix element.
    */
+  /// \{
   double& operator()(const size_t i, const size_t j);
-    
+  double& at(const size_t i, const size_t j);
+  /// \}
+
   /// Constant version of previous function.
+  /// \{
   const double& operator()(const size_t i, const size_t j) const;
+  const double& at(const size_t i, const size_t j) const;
+  /// \}
 
   /*! \brief Indexing for vectorisation.
    *
@@ -437,6 +443,7 @@
 
 //Inlined functions
 namespace kwantix{
+
 inline double& matrix::operator()(const size_t i, const size_t j){
   try{
     if(LUfactored)
@@ -456,8 +463,15 @@
     throw exc;
   }
 }
+
+inline double& matrix::at(const size_t i, const size_t j)
+{
+  return operator()(i,j);
+}
+
 inline const double& matrix::operator()(const size_t i,  
-                                        const size_t j) const{
+                                        const size_t j) const
+{
   try{
     return(*gsl_matrix_const_ptr(A,i-1,j-1));
   }
@@ -470,6 +484,11 @@
   }
 }
 
+inline const double& matrix::at(const size_t i, const size_t j) const
+{
+  return operator()(i,j);
+}
+
 inline double& vector::operator()(const size_t i) {
   try{
     return *gsl_vector_ptr(x,i-1);
@@ -491,5 +510,6 @@
     throw exc;
   }
 }
+
 }
 #endif //__LINALG_HPP__

--- a/src/interp_values.cpp
+++ b/src/interp_values.cpp
@@ -1,6 +1,4 @@
 #include "include/interp_values.hpp"
-//debug
-#include <iostream>
 
 namespace kwantix{
 
@@ -132,14 +130,13 @@
                  const map<point, kwantix::vector>& normals_in, 
                  size_t n_in) : rbfs_hash(rbfs_hash_in), n(n_in)
 {
-  map<point, kwantix::vector>::const_iterator I;
   size_t rows = normals_in.size();
   size_t cols = normals_in.begin() -> second.size();
   slice s(1,cols);
   size_t i = 1;
-  matrix N_temp(rows,cols);
-  for(I = normals_in.begin(); I != normals_in.end(); I++){
-    N_temp(i,s) = (I->second)/(norm(I->second));
+  shared_ptr<matrix> N_temp(new matrix(rows,cols));
+  for(auto I = normals_in.begin(); I != normals_in.end(); I++){
+    N_temp -> operator()(i,s) = (I->second)/(norm(I->second));
     i++;
   }
   N = N_temp;
@@ -147,8 +144,8 @@
 
 bdry_values normals::operator()(size_t k) const
 {
-  slice s(1,N.rows());
-  return bdry_values(N(s,k),rbfs_hash,n,m);
+  slice s(1,N -> rows());
+  return bdry_values(N -> operator()(s,k),rbfs_hash,n,m);
 }
 
 bool normals::operator==(const normals& M) const

--- a/src/interpolator.cpp
+++ b/src/interpolator.cpp
@@ -153,7 +153,7 @@
     exc.file = __FILE__;
     throw exc;
   }
-  coeffs = M.inv(values.v);
+  coeffs = M -> inv(values.v);
    
   //Precompute the relevant values
   recompute_values_vec();
@@ -192,7 +192,7 @@
   }
     
   //Now define the matrix to be inverted...
-  matrix Mtemp(n+m,n+m);
+  shared_ptr<matrix> Mtemp(new matrix(m+n,m+n));
   M = Mtemp;
   shared_ptr<const linear_diff_op2> L = thebvp -> get_linear_diff_op2(); 
 
@@ -201,7 +201,7 @@
   I = interior.begin();
   for(size_t i = 1; i <= n; i++){
     for(size_t j = 1; j <= n+m; j++)
-      M(i,j) = L -> at(rbfs[j-1], *I);
+      M -> at(i,j) = L -> at(rbfs[j-1], *I);
     I++;
   }
     
@@ -209,7 +209,7 @@
   J = normals.begin();
   for(size_t i = n+1; i <= n+m; i++){
     for(size_t j = 1; j <= n+m; j++)
-      M(i,j) = B -> at(rbfs[j-1], J->first, J->second);
+      M -> at(i,j) = B -> at(rbfs[j-1], J->first, J->second);
     J++;
   }
    
@@ -253,18 +253,18 @@
     std::vector<size_t> alpha; //empty vector
       
     //Precompute the RBFs at all points of the domain
-    matrix phis(n+m,n+m);
+    shared_ptr<matrix> phis(new matrix(n+m,n+m));
     set<point>::iterator I;
     size_t i;
     for(I = thebvp -> get_domain () -> get_interior().begin(), i = 1; 
         i <= n; i++, I++)
       for(size_t j = 0; j < n+m; j++)
-        phis(i,j+1) = rbfs[j].at(*I);
+        phis -> at(i,j+1) = rbfs[j].at(*I);
       
     for(I = thebvp -> get_domain() -> get_boundary().begin(), i=n+1; 
         i <= n+m; i++, I++)
       for(size_t j = 0; j < n+m; j++)
-        phis(i,j+1) = rbfs[j].at(*I);
+        phis -> at(i,j+1) = rbfs[j].at(*I);
       
     precomp_rbfs[alpha] = phis;
     computecoeffs();
@@ -287,19 +287,19 @@
     for(size_t k = 1; k <= thebvp -> get_domain() -> get_dimension(); 
         k++){
       std::vector<size_t> alpha(k); alpha[k-1]++;
-      matrix phis(n+m,n+m);
+      shared_ptr<matrix> phis(new matrix(n+m,n+m));
       {
         set<point>::iterator I;
         size_t i;
         for(I = thebvp -> get_domain() -> get_interior().begin(), 
               i = 1; i <= n; i++, I++)
           for(size_t j = 0; j < n+m; j++)
-            phis(i,j+1) = rbfs[j].d(*I,k);
+            phis -> at(i,j+1) = rbfs[j].d(*I,k);
       
         for(I = thebvp -> get_domain() -> get_boundary().begin(), 
               i=n+1; i <= n+m; i++, I++)
           for(size_t j = 0; j < n+m; j++)
-            phis(i,j+1) = rbfs[j].d(*I,k);
+            phis -> at(i,j+1) = rbfs[j].d(*I,k);
       }
 
       precomp_rbfs[alpha] = phis;
@@ -317,21 +317,24 @@
   if(!d2_precomp){
     size_t dim = thebvp -> get_domain() -> get_dimension();
     //precompute all second derivatives
-    for(size_t k1 = 1; k1 <= dim; k1++){
-      for(size_t k2 = k1; k2 <= dim; k2++){
+    for(size_t k1 = 1; k1 <= dim; k1++)
+    {
+      for(size_t k2 = k1; k2 <= dim; k2++)
+      {
         std::vector<size_t> alpha(k1>k2?k1:k2); alpha[k1-1]++; alpha[k2-1]++;
-        matrix phis(n+m,n+m);
+        shared_ptr<matrix> phis(new matrix(n+m,n+m));
+        
+        size_t i;
         set<point>::iterator I;
-        size_t i;
         for(I = thebvp -> get_domain () -> get_interior().begin(), i = 1; 
             i <= n; i++, I++)
           for(size_t j = 0; j < n+m; j++)
-            phis(i,j+1) = rbfs[j].d2(*I,k1,k2);
+            phis -> at(i,j+1) = rbfs[j].d2(*I,k1,k2);
       
         for(I = thebvp -> get_domain() -> get_boundary().begin(), i=n+1; 
             i <= n+m; i++, I++)
           for(size_t j = 0; j < n+m; j++)
-            phis(i,j+1) = rbfs[j].d2(*I,k1,k2);
+            phis -> at(i,j+1) = rbfs[j].d2(*I,k1,k2);
       
         precomp_rbfs[alpha] = phis;
         computecoeffs();	  
@@ -600,26 +603,31 @@
 
 template<typename RBF>
 void interpolator<RBF>::set_bdry_values(const bdry_values& b_new){
-  if(!initted){
+  if(!initted)
+  {
     throw not_initted(__LINE__, __FILE__);
   }
+  if(!precomputed)
+  {
+    throw not_precomputed(__LINE__, __FILE__);
+  }
+
   std::vector<size_t> alpha; //empty, corresponds to evaluation
-  if(precomp_values_vec.find(alpha) == precomp_values_vec.end())
-    at();
   kwantix::vector rhs = precomp_values_vec[alpha];
   slice s1(1,n), s2(n+1,n+m);
   rhs(s2) = b_new.v;
-  coeffs = precomp_rbfs[alpha].inv(rhs);
+  coeffs = precomp_rbfs[alpha] -> inv(rhs);
     
   recompute_values_vec();
 }
 
 template<typename RBF>
-void interpolator<RBF>::recompute_values_vec(){
-  typename 
-    map<std::vector<size_t>, matrix>::const_iterator I;
-  for(I = precomp_rbfs.begin(); I != precomp_rbfs.end(); I++)
-    precomp_values_vec[I -> first] = (I -> second)*coeffs;
+void interpolator<RBF>::recompute_values_vec()
+{
+  for(auto I = precomp_rbfs.begin(); I != precomp_rbfs.end(); I++)
+  {
+    precomp_values_vec[I -> first] = *(I -> second)*coeffs;
+  }
       
   //Pointwise map evaluations we will not compute again.
   precomp_values.clear();
@@ -662,15 +670,13 @@
       I++;
     }
   }
-  coeffs = M.inv(rhs);
+  coeffs = M -> inv(rhs);
     
   //Precompute the values for the derivatives and evaluations that
   //have already been computed.
   {				
-    typename map<std::vector<size_t>, matrix>::const_iterator I;
-      
-    for(I = precomp_rbfs.begin(); I != precomp_rbfs.end(); I++){
-      kwantix::vector vals = (I->second)*coeffs;
+    for(auto I = precomp_rbfs.begin(); I != precomp_rbfs.end(); I++){
+      kwantix::vector vals = *(I->second)*coeffs;
       precomp_values_vec[I -> first] = vals;
       precomp_values[I -> first] = valsvec2map(vals);
     }

--- a/src/linalg.cpp
+++ b/src/linalg.cpp
@@ -51,7 +51,7 @@
   if(SVDfactored){
     SVD = gsl_vector_alloc(M.SVD->size);
     gsl_vector_memcpy(SVD,M.SVD);
-  }       
+  } 
 }
 
 matrix matrix::operator=(const matrix& M){
@@ -114,7 +114,8 @@
   
 
   
-vector_view matrix::operator()(const size_t  i, const slice &b){
+vector_view matrix::operator()(const size_t  i, const slice &b)
+{
   vector_view x_sub(A,i,b);
 
   if(LUfactored)
@@ -126,12 +127,15 @@
   SVDfactored = false;
   return x_sub;
 }
-const vector_view matrix::operator()(const size_t  i, const slice &b) const{
+
+const vector_view matrix::operator()(const size_t  i, const slice &b) const
+{
   vector_view x_sub(A,i,b);
   return x_sub;
 }
   
-vector_view matrix::operator()(const slice &a, const size_t j){
+vector_view matrix::operator()(const slice &a, const size_t j)
+{
   vector_view x_sub(A,a,j);
 
   if(LUfactored)

--- a/src/main-sw-rk4.cpp
+++ b/src/main-sw-rk4.cpp
@@ -19,7 +19,7 @@
 
 using namespace kwantix;
 
-#define RBF_TYPE kwantix::multiquadric
+typedef kwantix::inverse_multiquadric RBF_TYPE;
 
 //Will solve the system
 //
@@ -160,12 +160,12 @@
     plotter.begin_interaction();
 
     //main loop, timestepping with RK4.
-    size_t maxiter = 1000;
+    size_t maxiter = 200;
     for(size_t i = 1; i <= maxiter; i++){ 
       cout << "Now on iteration #" << i << endl;
       plotter.update_values(h0.at());
       plotter.plot();
-      if(true)//i % 1 == 0)
+      if(false)//i % 1 == 0)
       {
         save_interp(Omega,u0,i,"u");
         save_interp(Omega,v0,i,"v");
@@ -294,8 +294,8 @@
   ofstream ofs(filename.c_str());    
   {
     slice s(1,2);
-    matrix M(Omega -> get_interior().size() +
-             Omega -> get_boundary().size(), 3);
+    static matrix M(Omega -> get_interior().size() +
+                    Omega -> get_boundary().size(), 3);
     size_t k = 1;
     for(set<point>::const_iterator i = Omega -> get_interior().begin();
         i != Omega -> get_interior().end(); i++){