diff --git a/Matrix.h b/Matrix.h
index ed6cc3f..8f9de2d 100644
--- a/Matrix.h
+++ b/Matrix.h
@@ -34,6 +34,210 @@ void print(const Slice& sl) {
   if (sl.ndim == 2) printf("Slice: i0 = %d  i1 = %d  j0 = %d  j1 = %d\n", sl.i0, sl.i1, sl.j0, sl.j1);
 }
 
+template <typename T>
+class VectorT {
+  // the dimensions
+  unsigned int _dim0;
+  // the data pointer
+  std::shared_ptr<T> _sp;
+  T* _p;
+  // is the matrix currently allocated?
+  bool _allocated;
+
+  friend VectorT copy(const VectorT& t, bool empty = false) {
+    VectorT r;
+    if (t._allocated) {
+      r._dim0 = t._dim0;
+      int sz = r.size();
+      r._p = new T[sz];
+      r._sp = std::shared_ptr<T>(r._p, [](T *p) {delete[] p;});
+      if (empty) {
+        for (int i = 0; i < sz; i++) {
+          r._p[i] = 0.0;
+        }
+      }
+      else {
+        for (int i = 0; i < sz; i++) {
+          r._p[i] = t._p[i];
+        }
+      }
+      r._allocated = t._allocated;
+    }
+    return r;
+  }
+
+public:
+  VectorT() 
+    :  _dim0(0), _p(0),_allocated(false) {}
+  
+  virtual ~VectorT() {
+    _allocated = false;
+    _dim0 = 0;
+  }
+  
+  VectorT& operator=(const VectorT& t) {
+    if (this != &t) {
+      _dim0 = t._dim0;
+      _p = t._p;
+      _sp = t._sp;
+      _allocated = t._allocated;
+    }
+    return *this;
+  }
+
+  VectorT(const VectorT& t) {
+    (*this) = t;
+  }
+
+  VectorT(const std::vector<T>& v) {
+    create(v.size());  
+    for (auto i = 0; i < v.size(); i++) _p[i] = v[i];
+  }
+
+  VectorT(int k) {
+    create(k);  
+    empty();
+  }
+
+  void create(int d0) {
+    // dims
+    _dim0 = d0;
+    // allocation
+    _p = new T[d0];
+    _sp = std::shared_ptr<T>(_p, [](T *p) {delete[] p;});
+    _allocated = true;
+  }
+
+  int size() const {
+    return _dim0;
+  }
+
+  // set all elements to zero
+  void empty() {
+    int sz = this->size();
+    for (int i = 0; i < sz; i++) {
+      _p[i] = T(0);
+    }
+  }
+
+  // set all elements to value
+  void value(T val) {
+    int sz = this->size();
+    for (int i = 0; i < sz; i++) {
+      _p[i] = val;
+    }
+  }
+
+  // pointer at a base position
+  T* ptr() {
+    return _p;
+  }
+
+  // pointer at a base position
+  const T* ptr() const {
+    return _p;
+  }
+
+  // pointer at a given index
+  T* ptr(int i0) {
+    return &_p[i0];
+  }
+
+  // pointer at a given index
+  const T* ptr(int i0) const {
+    return &_p[i0];
+  }
+
+  // index operator
+  T& operator[](int i) {
+    return _p[i];
+  }
+
+  // index operator
+  T& operator[](int i) const {
+    return _p[i];
+  }
+
+  // perform inplace a*(*this) + b*t
+  void gaxpy(const T& a, const VectorT& t, const T& b) {
+    int sz1 = this->size();
+    int sz2 = t.size();
+    assert(sz1 == sz2);
+    for (int i = 0; i < sz1; i++) _p[i] = a*_p[i]+b*t._p[i];
+  }
+
+  // perform out-of-place a*(*this) + b*t
+  VectorT gaxpy_oop(const T& a, const VectorT& t, const T& b) const {
+    int sz1 = this->size();
+    int sz2 = t.size();
+    assert(sz1 == sz2);
+    VectorT r = copy(*this, true);
+    for (int i = 0; i < sz1; i++) r._p[i] = a*_p[i]+b*t._p[i];
+    return r;
+  }
+
+  // addition
+  VectorT operator+(const VectorT& t) const {
+    return gaxpy_oop(1.0, t, 1.0);
+  }
+
+  // subtraction
+  VectorT operator-(const VectorT& t) const {
+    return gaxpy_oop(1.0, t, -1.0);
+  }
+
+  // inplace scaling by a constant 
+  void scale(T a) {
+    int sz = this->size();
+    for (int i = 0; i < sz; i++) _p[i] *= a;
+  }
+
+  // inner like a vector 
+  T inner(const VectorT& t) const {
+    int sz1 = this->size();
+    int sz2 = t.size();
+    assert(sz1 == sz2);
+    T rval = 0.0;
+    for (int i = 0; i < sz1; i++) rval += _p[i]*t._p[i]; 
+    return rval;
+  }
+
+  // two norm of the vector
+  T norm2() const {
+    int sz = this->size();
+    T rval = 0.0;
+    for (int i = 0; i < sz; i++) rval += _p[i]*_p[i]; 
+    return std::sqrt(rval);
+  }
+
+  // normalize the values according to the two-norm of the vector
+  void normalize() {
+    T s = this->norm2();
+    this->scale(1./s); 
+  }
+
+  VectorT<T> slice(int i0, int i1) {
+    VectorT<T> v;
+    int sz = i1-i0+1;
+    assert(sz > 0);
+    v.create(i1-i0+1);
+    for (int i = 0; i < sz; i++) 
+      v[i] = _p[i+i0];
+    return v;
+  }
+
+  T* begin() {
+    return &_p[0];
+  }
+
+  T* end() {
+    return &_p[size()-1];
+  }
+};
+
+typedef VectorT<double> real_vector;
+typedef VectorT<std::complex<double> > complex_vector;
+
 template <typename T>
 class MatrixT {
 private:
@@ -302,11 +506,25 @@ class MatrixT {
     return R;
   }
 
+  // multiply with a VectorT
+  VectorT<T> operator* (const VectorT<T>& v) const {
+    int vsz = v.size();
+    VectorT<T> r(vsz);
+    assert(vsz == _dim1);
+    for (int i = 0; i < _dim0; i++) {
+      T s = T(0);
+      for (int j = 0; j < _dim1; j++) {
+        s += _p[i*_dim0+j]*v[j];
+      }
+      r[i] = s;
+    }
+    return r;
+  }
+
   // multiply with a C++ std::vector
   std::vector<T> operator* (const std::vector<T>& v) const {
     int vsz = v.size();
     std::vector<T> r(vsz);
-    printf("vsz: %d   _dim1: %d\n", vsz, _dim0);
     assert(vsz == _dim1);
     for (int i = 0; i < _dim0; i++) {
       T s = T(0);
@@ -415,7 +633,7 @@ MatrixT<Q> constant(int d0, int d1, double val) {
 // create an constant matrix
 template <typename Q>
 MatrixT<Q> ones(int d0, int d1) {
-  return constant(d0, d1, T(1.0));
+  return constant(d0, d1, Q(1.0));
 }
 
 template <typename Q>
diff --git a/function1d.h b/function1d.h
index a644f83..8db92b9 100644
--- a/function1d.h
+++ b/function1d.h
@@ -3,7 +3,7 @@
 #include "twoscale.h"
 #include "gauss_legendre.h"
 
-using Vector = vector<double>;
+using Vector = real_vector;
 
 double normf(Vector v) {
   auto s = 0.0;
@@ -16,6 +16,7 @@ class Function1D {
   int k = 8;
   double thresh = 1e-6;
   int maxlevel = 30;
+  int initiallevel = 2;
   std::map<Key,Vector> tree;
   Vector quad_x;
   Vector quad_w;
@@ -28,11 +29,19 @@ class Function1D {
 
 public:
   Function1D(int k, double thresh, int maxlevel = 30) 
-   : k(k), thresh(thresh) {
+   : k(k), thresh(thresh), maxlevel(maxlevel) {
     init_twoscale(k);
     init_quadrature(k);
   }
 
+  Function1D(double (*f) (double), int k, double thresh, int maxlevel = 30, int initiallevel = 2) 
+   : k(k), thresh(thresh), maxlevel(maxlevel), initiallevel(initiallevel) {
+    init_twoscale(k);
+    init_quadrature(k);
+    int ntrans = std::pow(2, initiallevel);
+    for (auto l = 0; l < ntrans; l++) refine(f, initiallevel, l);
+  }
+
   ~Function1D() {}
 
   void init_twoscale(int k) {
@@ -44,6 +53,7 @@ class Function1D {
     quad_x = gauss_legendre_x(order);
     quad_w = gauss_legendre_w(order);
     quad_npts = quad_w.size();
+
     quad_phi  = zeros<double>(quad_npts, k);
     quad_phiT = zeros<double>(k, quad_npts);
     quad_phiw = zeros<double>(quad_npts, k);
@@ -52,7 +62,7 @@ class Function1D {
       for (auto m = 0; m < k; m++) {
         quad_phi(i,m) = p[m];
         quad_phiT(m,i) = p[m];
-        quad_phiw[i,m] = quad_w[i]*p[m];
+        quad_phiw(i,m) = quad_w[i]*p[m];
       }
     }
   }
@@ -72,19 +82,35 @@ class Function1D {
   }
 
   void refine(double (*f)(double), int n, int l) {
-    printf("refine n: %d     l: %d\n", n, l);
+    printf("\nrefine---> n: %d     l: %d\n", n, l);
     Vector s0 = project_box(f, n+1, 2*l);
+    printf("  computed s0\n");
     Vector s1 = project_box(f, n+1, 2*l+1);
+    printf("  computed s1\n");
     Vector s(2*k);
     for (auto i = 0; i < k; i++) s[i] = s0[i];
+    printf("  copied s0 to s\n");
     for (auto i = k; i < 2*k; i++) s[i+k] = s1[i];
+    printf("  copied s1 to s\n");
     Vector d = hg*s;
-    if (normf(Vector(d.begin()+k,d.begin()+2*k)) < thresh || n >= maxlevel-1) {
+    printf("  d = hg*s\n");
+    if (normf(Vector(d.slice(k,2*k-1))) < thresh || n >= maxlevel-1) {
       tree[Key(n+1,2*l)] = s0;
+      printf("set n+1 2*l coeff\n");
       tree[Key(n+1,2*l+1)] = s1;
+      printf("set n+1 2*l+1 coeff\n");
     } else {
       refine(f, n+1, 2*l);
       refine(f, n+1, 2*l+1);
     }
   }
+
+  void print_coeffs(int n, int l) {
+    auto s = tree[Key(n,l)];
+    printf("[%d, %d] (", n, l);
+    for (auto v : s) {
+      printf("%8.4f  ", v);
+    }
+    printf(")\n");
+  }
 };
diff --git a/gauss_legendre.h b/gauss_legendre.h
index 0ef5658..d4421db 100644
--- a/gauss_legendre.h
+++ b/gauss_legendre.h
@@ -1,7 +1,10 @@
 #ifndef GAUSS_LEGENDRE_H_
 #define GAUSS_LEGENDRE_H_
 
-#include <vector>
+#include "Matrix.h"
+
+//using GLVector = real_vector;
+using GLVector = std::vector<double>;
 
 #define GL_POINTS0 = {}
 #define GL_WEIGHTS0 = {}
@@ -38,53 +41,53 @@
 #define GL_POINTS11 = {0.98911432907302843, 0.94353129988404771, 0.86507600278702468, 0.75954806460340585, 0.63477157797617245, 0.50000000000000000, 0.36522842202382755, 0.24045193539659410, 0.13492399721297532, 0.05646870011595234, 0.01088567092697151}
 #define GL_WEIGHTS11 = { 0.02783428355808731, 0.06279018473245226, 0.09314510546386703, 0.11659688229599521, 0.13140227225512346, 0.13646254338895031, 0.13140227225512346, 0.11659688229599521, 0.09314510546386703, 0.06279018473245226, 0.02783428355808731}
 
-std::vector<double> gauss_legendre_x(int k) {
+GLVector gauss_legendre_x(int k) {
   if (k == 0) {
-    std::vector<double> x GL_POINTS0;
+    GLVector x GL_POINTS0;
     return x;
   }
   else if (k == 1) {
-    std::vector<double> x GL_POINTS1;
+    GLVector x GL_POINTS1;
     return x;
   }
   else if (k == 2) {
-    std::vector<double> x GL_POINTS2;
+    GLVector x GL_POINTS2;
     return x;
   }
   else if (k == 3) {
-    std::vector<double> x GL_POINTS3;
+    GLVector x GL_POINTS3;
     return x;
   }
   else if (k == 4) {
-    std::vector<double> x GL_POINTS4;
+    GLVector x GL_POINTS4;
     return x;
   }
   else if (k == 5) {
-    std::vector<double> x GL_POINTS5;
+    GLVector x GL_POINTS5;
     return x;
   }
   else if (k == 6) {
-    std::vector<double> x GL_POINTS6;
+    GLVector x GL_POINTS6;
     return x;
   }
   else if (k == 7) {
-    std::vector<double> x GL_POINTS7;
+    GLVector x GL_POINTS7;
     return x;
   }
   else if (k == 8) {
-    std::vector<double> x GL_POINTS8;
+    GLVector x GL_POINTS8;
     return x;
   }
   else if (k == 9) {
-    std::vector<double> x GL_POINTS9;
+    GLVector x GL_POINTS9;
     return x;
   }
   else if (k == 10) {
-    std::vector<double> x GL_POINTS10;
+    GLVector x GL_POINTS10;
     return x;
   }
   else if (k == 11) {
-    std::vector<double> x GL_POINTS11;
+    GLVector x GL_POINTS11;
     return x;
   }
   else {
@@ -92,53 +95,53 @@ std::vector<double> gauss_legendre_x(int k) {
   }
 }
 
-std::vector<double> gauss_legendre_w(int k) {
+GLVector gauss_legendre_w(int k) {
   if (k == 0) {
-    std::vector<double> w GL_WEIGHTS0;
+    GLVector w GL_WEIGHTS0;
     return w;
   }
   else if (k == 1) {
-    std::vector<double> w GL_WEIGHTS1;
+    GLVector w GL_WEIGHTS1;
     return w;
   }
   else if (k == 2) {
-    std::vector<double> w GL_WEIGHTS2;
+    GLVector w GL_WEIGHTS2;
     return w;
   }
   else if (k == 3) {
-    std::vector<double> w GL_WEIGHTS3;
+    GLVector w GL_WEIGHTS3;
     return w;
   }
   else if (k == 4) {
-    std::vector<double> w GL_WEIGHTS4;
+    GLVector w GL_WEIGHTS4;
     return w;
   }
   else if (k == 5) {
-    std::vector<double> w GL_WEIGHTS5;
+    GLVector w GL_WEIGHTS5;
     return w;
   }
   else if (k == 6) {
-    std::vector<double> w GL_WEIGHTS6;
+    GLVector w GL_WEIGHTS6;
     return w;
   }
   else if (k == 7) {
-    std::vector<double> w GL_WEIGHTS7;
+    GLVector w GL_WEIGHTS7;
     return w;
   }
   else if (k == 8) {
-    std::vector<double> w GL_WEIGHTS8;
+    GLVector w GL_WEIGHTS8;
     return w;
   }
   else if (k == 9) {
-    std::vector<double> w GL_WEIGHTS9;
+    GLVector w GL_WEIGHTS9;
     return w;
   }
   else if (k == 10) {
-    std::vector<double> w GL_WEIGHTS10;
+    GLVector w GL_WEIGHTS10;
     return w;
   }
   else if (k == 11) {
-    std::vector<double> w GL_WEIGHTS11;
+    GLVector w GL_WEIGHTS11;
     return w;
   }
   else {
diff --git a/test_function.cc b/test_function.cc
index f55ddd1..666850b 100644
--- a/test_function.cc
+++ b/test_function.cc
@@ -9,7 +9,7 @@ double func(double x) {
 }
 
 int main(int argc, char** argv) {
-  Function1D f(8, 1e-6, 10);
-  f.refine(func, 2, 1); 
+
+  Function1D f(func, 8, 1e-6, 3);
   return 0;
 }
diff --git a/tree.h b/tree.h
index 7f291f6..c1cda39 100644
--- a/tree.h
+++ b/tree.h
@@ -1,20 +1,42 @@
 #include <map>
+#include <cmath>
 #include "Matrix.h"
 
 struct Key {
   int n;
   int l;
+  Key() : n(0), l(0) {}
   Key(int n, int l) : n(n), l(l) {}
-  bool operator< (const Key& key) const {
-    if ((n < key.n) || (n == key.n && l < key.l)) return true;
-    return false;
+
+  bool operator< (const Key &k) const {
+    printf("n: %d  l: %d    k.n: %d  k.l: %d\n", n, l, k.n, k.l);
+    return ((1<<n)+l < (k.n<<n)+k.l);
   }
+
+  // bool operator< (const Key& key) const {
+  //   if ((n < key.n) || (n == key.n && l < key.l)) return true;
+  //   return false;
+  // }
 };
 
-class Tree {
-private:
-  std::map<Key,real_matrix> _map;
-public:
-  Tree() {}
-  ~Tree() {}
+struct MyHashCompare {
+    static size_t hash( const Key& x ) {
+        size_t prime = 31;
+        size_t result = 1;
+        result = prime * result + x.l;
+        result = prime * result + x.n;
+        return result;
+    }
+    //! True if strings are equal
+    static bool equal( const Key& x, const Key& y ) {
+        return ((x.n == y.n) && (x.l == y.l));
+    }
 };
+
+// class Tree {
+// private:
+//   std::map<Key,real_matrix,MyHashCompare> _map;
+// public:
+//   Tree() {}
+//   ~Tree() {}
+// };