diff --git a/@treefun2/balancef.m b/@treefun2/balancef.m
new file mode 100644
index 0000000..cbdd8f6
--- /dev/null
+++ b/@treefun2/balancef.m
@@ -0,0 +1,730 @@
+function g = balancef(f)
+% fortran vol_tree_fix_lr
+% assume square boxes? therefore scale is computed from dom(1:2)
+% nd needs to be fixed later...
+% flatNeighbors and leafNeighbors also need to be updated in the future
+%
+
+nd = 1; % 
+dpars = zeros(1000,1);
+ipars = zeros(100,1);
+zpars = zeros(10,1);
+
+if ( isempty(f) )
+    g = f;
+    return
+end
+
+% dada
+dom = f.domain(:,1)';
+scale = (dom(2)-dom(1))/2;
+norder = f.n;
+ndim = numel(dom(1:end/2));
+
+% convert to vol_tree_build output format
+iper = 0;
+npbox = norder^ndim;
+grid = zeros(ndim,norder^2);
+nbmax = 2^ndim*f.id(end); % is 2^ndim enough, since some box needs to be refined more than once...
+nlmax = f.height(1); % is this ok?
+centers0 = 1/2*(f.domain(1:2:end,:)+f.domain(2:2:end,:));
+centers = zeros(ndim,nbmax); centers(:,1:size(centers0,2)) = centers0;
+nlevels = f.height(1);
+nboxes = f.id(end);
+boxsize = 2*scale./2.^(0:nlmax); 
+laddr = zeros(2,nlmax+1); 
+for ilev=0:nlmax
+  laddr(1,ilev+1) = sum(f.level<ilev)+1;
+  laddr(2,ilev+1) = sum(f.level<=ilev);
+end
+ilevel = zeros(nbmax,1); ilevel(1:nboxes) = f.level;
+iparent = zeros(nbmax,1); iparent(1:nboxes) = f.parent;
+nchild = zeros(nbmax,1); nchild(1:nboxes) = sum(f.children>0);
+ichild = zeros(2^ndim,nbmax); ichild(:,1:nboxes) = f.children; ichild(ichild==0) = -1; % -1 convention in fortran tree
+
+% this does not belong to fortran tree, but needed for treefun plot
+coeffs = cell(1,nbmax); coeffs(1:nboxes) = f.coeffs(1:nboxes);
+
+% call computecoll to get vol_tree_fix_lr coll info
+nnbors0 = zeros(nboxes,1); nbors0 = zeros(3^ndim,nboxes);
+[nnbors0,nbors0] = ...
+  computecoll(ndim,nlevels,nboxes,laddr,boxsize,...
+                  centers(:,1:nboxes),iparent(1:nboxes),nchild(1:nboxes),ichild(:,1:nboxes),iper,...
+                  nnbors0,nbors0);
+nnbors = zeros(nbmax,1); nbors = zeros(3^ndim,nbmax);
+nnbors(1:nboxes) = nnbors0; nbors(:,1:nboxes) = nbors0;
+
+% now fix lr (with coeffs, which is different from)
+[centers,nboxes,laddr,ilevel,iparent,nchild,ichild,nnbors,nbors,coeffs] = ...
+      vol_tree_fix_lr_treefun(ndim,iper,nd,dpars,zpars,ipars,...
+                               norder,npbox,grid,centers,nlevels,nboxes,boxsize,...
+                               nbmax,nlmax,laddr,ilevel,iparent,nchild,ichild,nnbors,nbors,...
+                               coeffs);
+
+centers = centers(:,1:nboxes);
+ilevel = ilevel(1:nboxes);
+iparent = iparent(1:nboxes);
+nchild = nchild(1:nboxes);
+ichild = ichild(:,1:nboxes);
+nnbors = nnbors(1:nboxes);
+nbors = nbors(:,1:nboxes);
+coeffs = coeffs(1:nboxes);
+
+% convert back to treefun format, try to follow treefun naming convention
+g = f;
+scl                 = boxsize(ilevel+1); % x,y,z all the same
+g.domain            = zeros(2*ndim,nboxes);
+g.domain(1:2:end,:) = centers - 1/2*scl;
+g.domain(2:2:end,:) = centers + 1/2*scl;
+g.level             = ilevel(:)';
+g.id                = 1:nboxes;
+g.parent            = iparent(:)';
+g.children          = ichild;
+g.height            = zeros(1,nboxes);
+g.height(nchild>0)  = 1;
+for k = length(g.id):-1:1
+  if ( ~(g.height(k)==0) )
+    g.height(k) = 1 + max(g.height(g.children(:,k)));
+  end
+end
+% redo some stuff
+g.col               = uint64(zeros(1,nboxes));
+g.row               = uint64(zeros(1,nboxes));
+g.morton            = uint64(zeros(1,nboxes));
+mc = 2^ndim;
+for ilev = 0:nlevels-1
+  for ibox = laddr(1,ilev+1):laddr(2,ilev+1)
+    if nchild(ibox) > 0
+      if ndim == 2
+        j              = 1;
+        jbox           = g.children(j,ibox);
+        g.col(jbox)    = 2*g.col(ibox);
+        g.row(jbox)    = 2*g.row(ibox);
+        g.morton(jbox) = cartesian2morton(g.col(jbox), g.row(jbox));
+        j              = 2;
+        jbox           = g.children(j,ibox);
+        g.col(jbox)    = 2*g.col(ibox) + 1;
+        g.row(jbox)    = 2*g.row(ibox);
+        g.morton(jbox) = cartesian2morton(g.col(jbox), g.row(jbox));
+        j              = 3;
+        jbox           = g.children(j,ibox);
+        g.col(jbox)    = 2*g.col(ibox);
+        g.row(jbox)    = 2*g.row(ibox) + 1;
+        g.morton(jbox) = cartesian2morton(g.col(jbox), g.row(jbox));
+        j              = 4;
+        jbox           = g.children(j,ibox);
+        g.col(jbox)    = 2*g.col(ibox) + 1;
+        g.row(jbox)    = 2*g.row(ibox) + 1;
+        g.morton(jbox) = cartesian2morton(g.col(jbox), g.row(jbox));
+        % keyboard
+      end
+    end
+    % g.col(ibox) = 2*g.col(iparent(ibox));
+    % g.row(ibox) = 2*g.row(iparent(ibox));
+  end
+end
+g.coeffs = coeffs;
+end
+
+function [centers,nboxes,laddr,ilevel,iparent,nchild,ichild,nnbors,nbors,coeffs] = vol_tree_fix_lr_treefun(ndim,iperiod,nd,dpars,zpars,ipars, ...
+                                                                                    norder,npbox,grid,centers,nlevels,nboxes,boxsize, ...
+                                                                                    nbmax,nlmax,laddr,ilevel,iparent,nchild,ichild,nnbors,nbors, ...
+                                                                                    coeffs)
+
+laddrtail=zeros(2,nlmax+1);
+% boxsize(0:nlmax),laddr(2,0:nlmax),laddrtail(2,0:nlmax)
+idx1 = 1;
+%
+bs0=boxsize(0+idx1);
+%
+mc=2^ndim;
+mnbors=3^ndim;
+
+iflag = zeros(nbmax,1); iflag2 = zeros(nbmax,1);
+
+for i=1:nboxes
+   iflag(i) = 0;
+end
+
+% Flag boxes that violate level restriction by "1" Violatioin refers to any
+% box that is directly touching a box that is more than one level finer
+% Method:
+% 1) Carry out upward pass. For each box B, look at the colleagues of B's
+% grandparent 
+% 2) See if any of those colleagues are childless and in contact with B.
+%  
+% Note that we only need to get up to level two, as we will not find a
+% violation at level 0 and level 1 For such boxes, we set iflag(i) = 1
+% figure(1), clf,
+for ilev=nlevels:-1:2
+  distest = 1.05d0*(boxsize(ilev-1+idx1) + boxsize(ilev-2+idx1))/2.0d0;
+  for ibox = laddr(1,ilev+idx1):laddr(2,ilev+idx1) % look at current box
+    idad = iparent(ibox); % its parent box
+    igranddad = iparent(idad); % its grand parent box
+    
+    for i=1:nnbors(igranddad) % look at neighbor of grand parent
+       jbox = nbors(i,igranddad);
+       if((nchild(jbox)==0) && (iflag(jbox)==0)) % if neighbor of grand parent has no child, and has not been flagged
+          ict = 0; % count overlapping dimension between ibox & its grand parent's neighbor
+          for k=1:ndim
+             dis = centers(k,jbox) - centers(k,idad);
+             if (iperiod==1)
+                dp1=dis+bs0;
+                dm1=dis-bs0;
+                if (abs(dis).gt.abs(dp1)), dis=dp1; end
+                if (abs(dis).gt.abs(dm1)), dis=dm1; end
+             end
+             if(abs(dis)<=distest), ict = ict + 1; end
+          end
+          if(ict==ndim)
+            iflag(jbox) = 1;
+            % figure(1),
+            % plotibox(centers(:,ibox),boxsize(ilev-1+idx1)/2);
+            % plotibox(centers(:,jbox),4*boxsize(ilev-1+idx1)/2);
+            % axis equal
+          end
+       end
+    end
+  end
+end
+
+% Find all boxes that need to be given a flag+ A flag+ box will be denoted
+% by setting iflag(box) = 2 This refers to any box that is not already
+% flagged and is bigger than and is contacting a flagged box or another box
+% that has already been given a flag +. It is found by performing an upward
+% pass and looking at the flagged box's parents colleagues and a flag+
+% box's parents colleagues and seeing if they are childless and present the
+% case where a bigger box is contacting a flagged or flag+ box.
+% figure(2), clf,
+%% if secondary violator, 1st send it downstairs during Step 4, then wait to be processed in Step 5 or no need to process at all, seems to be later, from the "worst-case" example 
+for ilev = nlevels:-1:1
+  distest = 1.05d0*(boxsize(ilev+idx1) + boxsize(ilev-1+idx1))/2.0d0;
+  for ibox = laddr(1,ilev+idx1):laddr(2,ilev+idx1)
+    if((iflag(ibox)==1) || (iflag(ibox)==2)) % look at current box
+      idad = iparent(ibox); % its parent box
+      for i=1:nnbors(idad) % look at neighbor of parent
+        jbox = nbors(i,idad);
+        if((nchild(jbox)==0) && (iflag(jbox)==0)) % if neighbor of parent has no child, and has not been flagged
+          ict = 0; % count overlapping dimension between ibox & its parent's neighbor
+          for k=1:ndim
+            dis = centers(k,jbox) - centers(k,ibox);
+            if (iperiod==1)
+              dp1=dis+bs0;
+              dm1=dis-bs0;
+              if (abs(dis)>abs(dp1)), dis=dp1; end
+              if (abs(dis)>abs(dm1)), dis=dm1; end
+            end
+            if(abs(dis)<=distest), ict = ict + 1; end
+          end
+          if(ict==ndim)
+            iflag(jbox) = 2;
+            % figure(2),
+            % plotibox(centers(:,ibox),boxsize(ilev-1+idx1)/2);
+            % plotibox(centers(:,jbox),4*boxsize(ilev-1+idx1)/2);
+            % axis equal
+          end
+        end
+      end
+    end
+  end
+end
+
+% Subdivide all flag and flag+ boxes. Flag all the children of flagged
+% boxes as flag++. Flag++ boxes are denoted by setting iflag(box) = 3. The
+% flag++ boxes need to be checked later to see which of them need further
+% refinement. While creating new boxes, we will need to update all the tree
+% structures as well. Note that all the flagged boxes live between levels 1
+% and nlevels - 2. We process the boxes via a downward pass. We first
+% determine the number of boxes that are going to be subdivided at each
+% level and everything else accordingly
+for ilev = 0:nlevels
+   laddrtail(1,ilev+idx1) = 0;
+   laddrtail(2,ilev+idx1) = -1;
+end
+
+nboxes0 = nboxes;
+%% Frank's thesis, Appendix A, Step 4, divide all of the primary and secondary violators once. I think
+% figure out iflag = 3 case
+for ilev = 1:nlevels-2
+
+  laddrtail(1,ilev+1+idx1) = nboxes+1;
+
+  nbloc = laddr(2,ilev+idx1)-laddr(1,ilev+idx1)+1;
+  
+  [iflag,centers,nboxes,ilevel,iparent,nchild,ichild,coeffs]  = vol_tree_refine_boxes_flag_treefun(ndim,iflag,nd,npbox,...
+                                   dpars,zpars,ipars,grid,nbmax,laddr(1,ilev+idx1),nbloc, ...
+                                   centers,boxsize(ilev+1+idx1),nboxes,ilev,ilevel,iparent,nchild,ichild, ...
+                                   coeffs);
+  
+  laddrtail(2,ilev+1+idx1) = nboxes;
+
+end
+
+ishuffle1 = zeros(nboxes,1);
+centerstmp = centers(:,1:nboxes); ileveltmp = ilevel(1:nboxes); iparenttmp = iparent(1:nboxes); nchildtmp = nchild(1:nboxes); ichildtmp = ichild(:,1:nboxes); iflagtmp = iflag(1:nboxes);
+coeffstmp = coeffs(1:nboxes);
+[centerstmp,laddr,ileveltmp,ishuffle1,iparenttmp,nchildtmp,ichildtmp,coeffstmp,iflagtmp] = ... 
+              vol_tree_reorg_treefun(ndim,nboxes,nd,npbox,centerstmp,nlevels,laddr,...
+                                      laddrtail,ileveltmp,ishuffle1,iparenttmp,nchildtmp,ichildtmp,coeffstmp,iflagtmp);
+centers(:,1:nboxes) = centerstmp; ilevel(1:nboxes) = ileveltmp; iparent(1:nboxes) = iparenttmp; nchild(1:nboxes) = nchildtmp; ichild(:,1:nboxes) = ichildtmp; iflag(1:nboxes) = iflagtmp;
+coeffs(1:nboxes) = coeffstmp;
+
+nnborstmp = nnbors(1:nboxes); nborstmp = nbors(:,1:nboxes);
+[nnborstmp,nborstmp] = computecoll(ndim,nlevels,nboxes,laddr,boxsize,...
+                                      centers(:,1:nboxes),iparent(1:nboxes),nchild(1:nboxes),...
+                                      ichild(:,1:nboxes),iperiod,nnborstmp,nborstmp);
+nnbors(1:nboxes) = nnborstmp; nbors(:,1:nboxes) = nborstmp;
+
+
+% Processing of flag and flag+ boxes is done Start processing flag++ boxes.
+% We will use a similar strategy as before. We keep checking the flag++
+% boxes that require subdivision if they still violate the level
+% restriction criterion, create the new boxes, append them to the end of
+% the list to begin with and in the end reorganize the tree structure. We
+% shall accomplish this via a downward pass as new boxes that get added in
+% the downward pass will also be processed simultaneously. We shall
+% additionally also need to keep on updating the colleague information as
+% we proceed in the downward pass
+% 
+% Reset the flags array to remove all the flag and flag+ cases. This is to
+% ensure reusability of the subdivide_flag routine to handle the flag++ case
+
+for ibox=1:nboxes
+  if(iflag(ibox)~=3), iflag(ibox) = 0; end
+end
+
+for ilev = 0:nlevels
+  laddrtail(1,ilev+idx1) = 0;
+  laddrtail(2,ilev+idx1) = -1;
+end
+
+% boxsize(0:nlmax),laddr(2,0:nlmax),laddrtail(2,0:nlmax)
+%% Frank's thesis, Appendix A, Step 5, at each level, for a descendant of a primary violator, test whether it is still in violation
+for ilev = 2:nlevels-2
+
+  % Step 1
+  iflagtmp = iflag(1:nboxes);
+  iflagtmp = vol_updateflags(ndim,iperiod,ilev,nboxes,nlevels, ...
+            laddr,nchild(1:nboxes),ichild(:,1:nboxes),nnbors(1:nboxes),nbors(:,1:nboxes),centers(:,1:nboxes),boxsize,iflagtmp);
+  iflag(1:nboxes) = iflagtmp;
+
+  iflagtmp = iflag(1:nboxes);
+  iflagtmp = vol_updateflags(ndim,iperiod,ilev,nboxes,nlevels, ...
+            laddrtail,nchild(1:nboxes),ichild(:,1:nboxes),nnbors(1:nboxes),nbors(:,1:nboxes),centers(:,1:nboxes),boxsize,iflagtmp);
+  iflag(1:nboxes) = iflagtmp;
+
+  % Step 2
+  laddrtail(1,ilev+1+idx1) = nboxes + 1;
+
+  nbloc = laddr(2,ilev+idx1)-laddr(1,ilev+idx1)+1;
+  [iflag,centers,nboxes,ilevel,iparent,nchild,ichild,coeffs] = vol_tree_refine_boxes_flag_treefun(ndim,iflag,nd,npbox, ...
+                                                       dpars,zpars,ipars,grid,nbmax,laddr(1,ilev+idx1),nbloc, ...
+                                                       centers,boxsize(ilev+1+idx1),nboxes,ilev,ilevel,iparent,nchild,ichild, ...
+                                                       coeffs);
+  
+  nbloc = laddrtail(2,ilev+idx1)-laddrtail(1,ilev+idx1)+1;
+  [iflag,centers,nboxes,ilevel,iparent,nchild,ichild,coeffs] = vol_tree_refine_boxes_flag_treefun(ndim,iflag,nd,npbox, ...
+                                                       dpars,zpars,ipars,grid,nbmax,laddrtail(1,ilev+idx1),nbloc, ...
+                                                       centers,boxsize(ilev+1+idx1),nboxes,ilev,ilevel,iparent,nchild,ichild, ...
+                                                       coeffs);
+  laddrtail(2,ilev+1+idx1) = nboxes;      
+
+  % Step 3
+  for ibox = laddrtail(1,ilev+1+idx1):laddrtail(2,ilev+1+idx1)
+    nnbors(ibox) = 0;
+    idad = iparent(ibox);
+    for i=1:nnbors(idad)
+      jbox = nbors(i,idad);
+      for j=1:mc
+        kbox = ichild(j,jbox);
+        if(kbox>0) 
+          ifnbor=1;
+          for k=1:ndim
+            dis=centers(k,kbox)-centers(k,ibox);
+            if (iperiod==1) 
+              dp1=dis+bs0;
+              dm1=dis-bs0;
+              if (abs(dis)>abs(dp1)), dis=dp1; end
+              if (abs(dis)>abs(dm1)), dis=dm1; end
+            end
+            if((abs(dis)>1.05*boxsize(ilev+1+idx1)))
+              ifnbor=0;
+              break
+            end
+          end
+          if (ifnbor==1)
+            nnbors(ibox) = nnbors(ibox)+1;
+            nbors(nnbors(ibox),ibox) = kbox;
+          end
+        end
+      end
+    end
+    % End of computing colleagues of box i
+  end
+  
+end
+
+ishuffle2 = zeros(nboxes,1);
+centerstmp = centers(:,1:nboxes); ileveltmp = ilevel(1:nboxes); iparenttmp = iparent(1:nboxes); nchildtmp = nchild(1:nboxes); ichildtmp = ichild(:,1:nboxes); iflagtmp = iflag(1:nboxes);
+coeffstmp = coeffs(1:nboxes);
+[centerstmp,laddr,ileveltmp,ishuffle2,iparenttmp,nchildtmp,ichildtmp,coeffstmp,iflagtmp] = ... 
+              vol_tree_reorg_treefun(ndim,nboxes,nd,npbox,centerstmp,nlevels,laddr,...
+                                      laddrtail,ileveltmp,ishuffle2,iparenttmp,nchildtmp,ichildtmp,coeffstmp,iflagtmp);
+centers(:,1:nboxes) = centerstmp; ilevel(1:nboxes) = ileveltmp; iparent(1:nboxes) = iparenttmp; nchild(1:nboxes) = nchildtmp; ichild(:,1:nboxes) = ichildtmp; iflag(1:nboxes) = iflagtmp;
+coeffs(1:nboxes) = coeffstmp;
+
+nnborstmp = nnbors(1:nboxes); nborstmp = nbors(:,1:nboxes);
+[nnborstmp,nborstmp] = computecoll(ndim,nlevels,nboxes,laddr,boxsize,...
+                                      centers(:,1:nboxes),iparent(1:nboxes),nchild(1:nboxes),...
+                                      ichild(:,1:nboxes),iperiod,nnborstmp,nborstmp);
+nnbors(1:nboxes) = nnborstmp; nbors(:,1:nboxes) = nborstmp;
+
+end
+
+function [iflag,centers,nbctr,ilevel,iparent,nchild,ichild,coeffs] = vol_tree_refine_boxes_flag_treefun(ndim,iflag,nd,npbox,dpars,zpars,ipars,grid,nboxes,ifirstbox,nbloc,centers,bs,nbctr,nlctr,ilevel,iparent,nchild,ichild,coeffs)
+
+% myfun = @(x,y) sin(x).*cos(y)+exp(x.*y);
+% myvals = myfun(2*(xx0-1/2),2*(yy0-1/2));
+% mycoeffs = treefun2.vals2coeffs(myvals);
+% mycvals1 = reshape(treefun2.coeffs2checkvals(mycoeffs,cxx1(:),cyy1(:)),[norder norder]);
+% mycvals1_2 = myfun(cxx1,cyy1);
+% mycvals2 = reshape(treefun2.coeffs2checkvals(mycoeffs,cxx2(:),cyy2(:)),[norder norder]);
+% mycvals3 = reshape(treefun2.coeffs2checkvals(mycoeffs,cxx3(:),cyy3(:)),[norder norder]);
+% mycvals4 = reshape(treefun2.coeffs2checkvals(mycoeffs,cxx4(:),cyy4(:)),[norder norder]);
+% figure(1),clf,surf(cxx1,cyy1,mycvals1),
+% hold on, surf(cxx2,cyy2,mycvals2)
+% hold on, surf(cxx3,cyy3,mycvals3)
+% hold on, surf(cxx4,cyy4,mycvals4)
+
+% additional
+persistent xx0 yy0 xxx0 yyy0 nstored
+norder=(npbox)^(1/ndim);
+n=norder;
+if ndim == 2
+  nalias = n;
+  nrefpts = n; % Sample at equispaced points to test error
+
+  if ( isempty(xx0) || isempty(xxx0) || n ~= nstored )
+      nstored = n;
+      [xx0, yy0] = chebpts2(nalias, nalias, [0 1 0 1]);
+      [xxx0, yyy0] = meshgrid(linspace(0, 1, nrefpts));
+  end
+  sclx = 2*bs;
+  scly = 2*bs;
+end
+
+if ndim == 3
+end
+
+%
+isgn=zeros(ndim,2^ndim);
+isgn=get_child_box_sign(ndim,isgn);
+
+mc=2^ndim;
+
+ilastbox = ifirstbox+nbloc-1;
+
+bsh = bs/2;
+
+isum=zeros(nbloc,1);
+
+isum = cumsum(iflag(ifirstbox:(ifirstbox+nbloc-1))>0); % there are iflag = 2 entries
+
+for i = 1:nbloc
+  ibox = ifirstbox + i-1;
+  if(iflag(ibox)>0)
+    nbl = nbctr + (isum(i)-1)*mc;
+    nchild(ibox) = mc;
+    % additional
+    if ndim == 2
+      xx = sclx*xx0 + (centers(1,ibox)-1/2*sclx);
+      yy = scly*yy0 + (centers(2,ibox)-1/2*scly);
+      vals = treefun2.coeffs2vals(coeffs{ibox});
+      % split into four chilx boxes
+      dom = [-1 1 -1 1];
+      xmid = mean(dom(1:2));
+      ymid = mean(dom(3:4));
+      j=1;
+      jbox = nbl+j;
+      cdom1              = [dom(1) xmid dom(3) ymid];
+      csclx1             = diff(cdom1(1:2));
+      cscly1             = diff(cdom1(3:4));
+      cxx1               = csclx1*xx0 + cdom1(1); 
+      cyy1               = cscly1*yy0 + cdom1(3); 
+      cvals1             = reshape(treefun2.coeffs2checkvals(coeffs{ibox},cxx1(:),cyy1(:)),[norder norder]);
+      ccoeffs1           = treefun2.vals2coeffs(cvals1);
+      coeffs{jbox}       = ccoeffs1;
+      j=2;
+      jbox = nbl+j;
+      cdom2              = [xmid dom(2) dom(3) ymid];
+      csclx2             = diff(cdom2(1:2));
+      cscly2             = diff(cdom2(3:4));
+      cxx2               = csclx2*xx0 + cdom2(1); 
+      cyy2               = cscly2*yy0 + cdom2(3); 
+      cvals2             = reshape(treefun2.coeffs2checkvals(coeffs{ibox},cxx2(:),cyy2(:)),[norder norder]);
+      ccoeffs2           = treefun2.vals2coeffs(cvals2);
+      coeffs{jbox}       = ccoeffs2;
+      j=3;
+      jbox = nbl+j;
+      cdom3              = [dom(1) xmid ymid dom(4)];
+      csclx3             = diff(cdom3(1:2));
+      cscly3             = diff(cdom3(3:4));
+      cxx3               = csclx3*xx0 + cdom3(1); 
+      cyy3               = cscly3*yy0 + cdom3(3); 
+      cvals3             = reshape(treefun2.coeffs2checkvals(coeffs{ibox},cxx3(:),cyy3(:)),[norder norder]);
+      ccoeffs3           = treefun2.vals2coeffs(cvals3);
+      coeffs{jbox}       = ccoeffs3;
+      j=4;
+      jbox = nbl+j;
+      cdom4              = [xmid dom(2) ymid dom(4)];
+      csclx4             = diff(cdom4(1:2));
+      cscly4             = diff(cdom4(3:4));
+      cxx4               = csclx4*xx0 + cdom4(1); 
+      cyy4               = cscly4*yy0 + cdom4(3); 
+      cvals4             = reshape(treefun2.coeffs2checkvals(coeffs{ibox},cxx4(:),cyy4(:)),[norder norder]);
+      ccoeffs4           = treefun2.vals2coeffs(cvals4);
+      coeffs{jbox}       = ccoeffs4;
+      %
+      coeffs{ibox}       = [];
+    end
+    for j=1:mc
+      jbox = nbl+j;
+      for k=1:ndim
+        centers(k,jbox) = centers(k,ibox)+isgn(k,j)*bsh;
+      end
+      iparent(jbox) = ibox;
+      nchild(jbox) = 0;
+      for l=1:mc
+        ichild(l,jbox) = -1;
+      end
+      ichild(j,ibox) = jbox;
+      ilevel(jbox) = nlctr+1; 
+      if(iflag(ibox)==1), iflag(jbox) = 3; end
+      if(iflag(ibox)==2), iflag(jbox) = 0; end
+    end
+  end
+end
+
+if(nbloc>0), nbctr = nbctr + isum(nbloc)*mc; end
+
+end
+
+function iflag = vol_updateflags(ndim,iperiod,curlev,nboxes,nlevels,laddr,nchild,ichild,nnbors,nbors,centers,boxsize,iflag)
+%
+%
+
+idx1 = 1;
+bs0=boxsize(0+idx1);
+
+mc=2^ndim;
+distest = 1.05d0*(boxsize(curlev+idx1) + boxsize(curlev+1+idx1))/2.0d0;
+
+for ibox = laddr(1,curlev+idx1):laddr(2,curlev+idx1)
+  if(iflag(ibox)==3)
+    iflag(ibox) = 0;
+    for i=1:nnbors(ibox)
+      jbox = nbors(i,ibox);
+      
+      for j=1:mc
+        kbox = ichild(j,jbox);
+        if(kbox>0)
+          if(nchild(kbox)>0)
+            ict = 0;
+            for k=1:ndim
+              dis = centers(k,kbox) - centers(k,ibox);
+              if (iperiod==1) 
+                dp1=dis+bs0;
+                dm1=dis-bs0;
+                if (abs(dis)>abs(dp1)), dis=dp1; end
+                if (abs(dis)>abs(dm1)), dis=dm1; end
+              end
+              if(abs(dis)<=distest), ict = ict + 1; end
+            end
+            if(ict==ndim)
+               iflag(ibox) = 1;
+               break;
+            end
+          end
+        end
+      end
+      if iflag(ibox) == 1
+       break;
+      end
+    end
+    % 1111       continue   
+  end
+end
+
+end
+
+function [nnbors,nbors] = computecoll(ndim,nlevels,nboxes,laddr,boxsize,centers,iparent,nchild,ichild,iperiod,nnbors,nbors)
+
+% laddr(2,0:nlevels)
+% boxsize(0:nlevels)
+
+idx1 = 1;
+
+bs0=boxsize(0+idx1);
+
+mc= 2^ndim;
+mnbors=3^ndim;
+
+for i=1:nboxes
+  nnbors(i) = 0;
+  for j=1:mnbors
+    nbors(j,i) = -1;
+  end
+end
+
+nnbors(1) = 1;
+nbors(1,1) = 1;
+for ilev = 1:nlevels
+  ifirstbox = laddr(1,ilev+idx1);
+  ilastbox = laddr(2,ilev+idx1);
+  
+  for ibox = ifirstbox:ilastbox
+    dad = iparent(ibox);
+    for i=1:nnbors(dad)
+      jbox = nbors(i,dad);
+      for j=1:mc
+        kbox = ichild(j,jbox);
+        if(kbox>0)
+          ifnbor=1;
+          for k=1:ndim
+            dis=abs(centers(k,kbox)-centers(k,ibox));
+            if (iperiod==1)
+              dp1=bs0-dis;
+              if (dp1<dis), dis=dp1; end
+            end
+            if (dis>1.05*boxsize(ilev+idx1))
+              ifnbor=0;
+              break
+            end
+          end
+             
+          if(ifnbor==1)
+            nnbors(ibox) = nnbors(ibox)+1;
+            nbors(nnbors(ibox),ibox) = kbox;
+          end
+        end
+      end
+    end
+  end
+end
+
+end
+
+function [centers,laddr,ilevel,ishuffle,iparent,nchild,ichild,coeffs,iflag] = vol_tree_reorg_treefun(ndim,nboxes,nd,npbox,centers,nlevels,laddr,laddrtail,ilevel,ishuffle,iparent,nchild,ichild,coeffs,iflag)
+%
+%
+
+% laddrtail(2,0:nlevels)
+% laddr(2,0:nlevels)
+% tladdr(2,0:nlevels)
+
+tladdr = zeros(2,nlevels+1);
+
+idx1 = 1;
+
+mc=2^ndim;
+
+tilevel=zeros(nboxes,1); tiparent=zeros(nboxes,1); tnchild=zeros(nboxes,1);
+tichild=zeros(mc,nboxes); tiflag=zeros(nboxes,1); iboxtocurbox=zeros(nboxes,1);
+tcenters=zeros(ndim,nboxes);
+
+for ilev = 0:nlevels
+  tladdr(1,ilev+idx1) = laddr(1,ilev+idx1);
+  tladdr(2,ilev+idx1) = laddr(2,ilev+idx1);
+end
+
+tcenters = centers;
+tilevel = ilevel;
+tiparent = iparent;
+tnchild = nchild;
+tichild = ichild;
+tcoeffs = coeffs;
+
+for ibox=1:nboxes
+  tiflag(ibox) = iflag(ibox);
+end
+
+for ilev = 0:1
+  for ibox = laddr(1,ilev+idx1):laddr(2,ilev+idx1)
+    iboxtocurbox(ibox) = ibox;
+  end
+end
+
+ilevptr=zeros(nlevels+1,1); ilevptr2=zeros(nlevels,1);
+
+ilevptr(2) = laddr(1,2+idx1);
+
+for ilev=2:nlevels
+  nblev = laddr(2,ilev+idx1)-laddr(1,ilev+idx1)+1; % original num of boxes at ilev
+  ilevptr2(ilev) = ilevptr(ilev) + nblev;
+  nblev = laddrtail(2,ilev+idx1)-laddrtail(1,ilev+idx1)+1; % additional number of boxes at ilev due to refinement
+  ilevptr(ilev+1) = ilevptr2(ilev) + nblev;
+end
+
+curbox = laddr(1,2+idx1);
+ishuffle = zeros(nboxes,1);
+for ilev=2:nlevels
+  laddr(1,ilev+idx1) = curbox;
+  for ibox = tladdr(1,ilev+idx1):tladdr(2,ilev+idx1)
+    ilevel(curbox) = tilevel(ibox);
+    nchild(curbox) = tnchild(ibox);
+    ishuffle(ibox) = curbox; % ibox info copied to curbox
+    for k=1:ndim
+      centers(k,curbox) = tcenters(k,ibox);
+    end
+    coeffs{curbox} = tcoeffs{ibox};
+    iflag(curbox) = tiflag(ibox);
+    iboxtocurbox(ibox) = curbox;
+    
+    curbox = curbox + 1;
+  end
+  for ibox = laddrtail(1,ilev+idx1):laddrtail(2,ilev+idx1)
+    ilevel(curbox) = tilevel(ibox);
+    ishuffle(ibox) = curbox; % ibox info copied to curbox
+    for k=1:ndim
+      centers(k,curbox) = tcenters(k,ibox);
+    end
+    coeffs{curbox} = tcoeffs{ibox};
+    nchild(curbox) = tnchild(ibox);
+    iflag(curbox) = tiflag(ibox);
+    iboxtocurbox(ibox) = curbox;
+    
+    curbox = curbox + 1;
+  end
+  laddr(2,ilev+idx1) = curbox-1;
+end
+
+for ibox=1:nboxes
+  if(tiparent(ibox)==-1), iparent(iboxtocurbox(ibox)) = -1; end
+  if(tiparent(ibox)>0)
+    iparent(iboxtocurbox(ibox)) = iboxtocurbox(tiparent(ibox)); end
+  for i=1:mc
+    if(tichild(i,ibox)==-1), ichild(i,iboxtocurbox(ibox)) = -1; end
+    if(tichild(i,ibox)>0)
+      ichild(i,iboxtocurbox(ibox)) = iboxtocurbox(tichild(i,ibox)); end
+  end
+end
+
+end
+
+function isgn = get_child_box_sign(ndim,isgn)
+
+mc = 2^ndim;
+for j=1:ndim
+  isgn(j,1)=-1;
+end
+
+for j=1:ndim
+  for i=1:2^(j-1):mc
+    if (i>1), isgn(j,i)=-isgn(j,i-2^(j-1)); end
+    for k=1:2^(j-1)-1
+      isgn(j,i+k)=isgn(j,i);
+    end
+  end
+end
+
+end
\ No newline at end of file
diff --git a/@treefun2/coeffs2checkvals.m b/@treefun2/coeffs2checkvals.m
new file mode 100644
index 0000000..536b2b4
--- /dev/null
+++ b/@treefun2/coeffs2checkvals.m
@@ -0,0 +1,29 @@
+function vals = coeffs2checkvals(coeffs,x,y)
+%
+%
+
+persistent Evalx Evaly
+[p,~,nd] = size(coeffs);
+ncheckpts = numel(x); % hopefully not too large...
+
+Evalx = ones(ncheckpts, p);
+Evaly = ones(ncheckpts, p);
+Evalx(:,2) = x(:);
+Evaly(:,2) = y(:);
+for k=3:p
+  Evalx(:,k) = 2*x(:).*Evalx(:,k-1)-Evalx(:,k-2);
+  Evaly(:,k) = 2*y(:).*Evaly(:,k-1)-Evaly(:,k-2);
+end
+% coeffs to value map
+% vals = zeros(nd,ncheckpts);
+% for k=1:ncheckpts
+%   tmp1 = permute(tensorprod(Evalx(k,:),coeffs,2,1),[2 1 3]);
+%   vals(:,k) = squeeze(permute(tensorprod(Evaly(k,:),tmp1,2,1),[2 1 3]));
+% end
+vals = zeros(nd,ncheckpts);
+for k=1:ncheckpts
+  tmp1 = permute(tensorprod(Evaly(k,:),coeffs,2,1),[2 1 3]);
+  vals(:,k) = squeeze(permute(tensorprod(Evalx(k,:),tmp1,2,1),[2 1 3]));
+end
+
+end
\ No newline at end of file
diff --git a/@treefun2/treefun2.m b/@treefun2/treefun2.m
index f3afa13..b6b6380 100644
--- a/@treefun2/treefun2.m
+++ b/@treefun2/treefun2.m
@@ -158,6 +158,7 @@
             % Now do level restriction
             if ( opts.balance )
                 f = balance(f);
+                % f = balancef(f);
             else
                 % Do a cumulative sum in reverse to correct the heights
                 for k = length(f.id):-1:1
@@ -254,6 +255,7 @@
         vals = coeffs2vals(coeffs);
         vals = coeffs2refvals(coeffs);
         refvals = chebvals2refvals(chebvals);
+        checkvals = coeffs2checkvals(coeffs,x,y);
 
     end
 
diff --git a/@treefun3/balance.m b/@treefun3/balance.m
new file mode 100644
index 0000000..9967a1c
--- /dev/null
+++ b/@treefun3/balance.m
@@ -0,0 +1,308 @@
+function g = balance(f)
+% based on @treefun2.balance
+% initial attempt in BBBBBBoston
+
+if ( isempty(f) )
+    g = f;
+    return
+end
+
+% Note: We assume IDs are stored in level order
+
+% Unpack data because MATLAB seems to be slow for class member access...
+fmorton  = f.morton(:);
+fcol     = f.col;
+frow     = f.row;
+fdep     = f.dep;
+fxy      = [fcol(:) frow(:) fdep(:)];
+flevel   = f.level;
+fheight  = f.height;
+fcoeffs  = f.coeffs;
+fnlevels = max(fheight) + 1;
+fnboxes  = length(f.id);
+
+% Generate level indices for f
+flevelIdx = zeros(fnlevels+1, 1);
+idx = 1;
+flevelIdx(idx) = 1;
+currentLevel = 0;
+for k = 2:fnboxes
+    if ( currentLevel ~= flevel(k) )
+        idx = idx+1;
+        flevelIdx(idx) = k;
+        currentLevel = flevel(k);
+    end
+end
+flevelIdx(end) = fnboxes+1;
+
+% Generate a set of 2:1-balanced Morton IDs level by level from the bottom
+% (leaves) to the top (root)
+S = [];
+T = cell(fnlevels, 1);
+for l = fnlevels-1:-1:1
+    idx = flevelIdx(l):flevelIdx(l+1)-1;
+    fmortonl = fmorton(idx);
+    N = unique([S ; fmortonl(fheight(idx) > 0)]);
+    S = mortonparent(mortoncolleagues(N, l));
+    [child1, child2, child3, child4, child5, child6, child7, child8] = mortonchildren(N);
+    T{l+1} = unique([child1 ; child2 ; child3 ; child4 ; ...
+                     child5 ; child6 ; child7 ; child8]);
+
+    % fxyl = fxy(idx,:);
+    % N = unique([S ; fxyl(fheight(idx)>0, :)], 'rows', 'stable');
+    % S = xyparent(xycolleagues(N, l));
+    % [childxy1, childxy2, childxy3, childxy4] = xychildren(N);
+    % %childxy = [childxy1 childxy2 childxy3 childxy4].';
+    % cxy = zeros(4*size(childxy1,1), 2, 'uint64');
+    % cxy(1:4:end,:) = childxy1;
+    % cxy(2:4:end,:) = childxy2;
+    % cxy(3:4:end,:) = childxy3;
+    % cxy(4:4:end,:) = childxy4;
+    % T{l+1} = unique(cxy, 'rows', 'stable');
+    % T{l+1} = mortonsort(T{l+1});
+end
+T{1} = uint64(0);
+%T{1} = uint64([0 0]);
+
+% Generate level indices for g
+gnlevels = fnlevels;
+gnboxes = 0;
+glevelIdx = zeros(gnlevels+1, 1);
+glevelIdx(1) = 1;
+for l = 1:gnlevels
+    gnboxes = gnboxes + size(T{l}, 1);
+    glevelIdx(l+1) = gnboxes + 1;
+end
+
+% First, make the tree from the Morton IDs
+gmorton = cell2mat(T).';
+[gcol, grow, gdep] = morton2cartesian(gmorton);
+%gmorton = [];
+%xy = cell2mat(T);
+%gcol = xy(:,1).';
+%grow = xy(:,2).';
+glevel    = zeros(1, gnboxes);
+gheight   = zeros(1, gnboxes);
+gparent   = zeros(1, gnboxes);
+gchildren = zeros(8, gnboxes);
+gcoeffs   = cell(1, gnboxes);
+
+for l = gnlevels:-1:2
+    id = glevelIdx(l);
+    jparent = 1;
+    ps = mortonparent(T{l});
+    %pxys = xyparent(T{l});
+    for k = 1:8:size(T{l},1)
+        ids = id:id+7;
+        glevel(ids) = l-1;
+        p = ps(k);
+        %pxy = pxys(k,:);
+        while ( T{l-1}(jparent) ~= p && jparent <= length(T{l-1}) )
+        %while ( ~all(T{l-1}(jparent,:) == pxy) && jparent <= size(T{l-1}, 1) )
+            jparent = jparent + 1;
+        end
+        if ( jparent > size(T{l-1}, 1) )
+            error('Couldn''t find parent node.');
+        end
+        pid = glevelIdx(l-1) + jparent - 1;
+        gparent(ids) = pid;
+        gchildren(:,pid) = ids;
+        gheight(pid) = 1 + max(gheight(ids));
+        id = id + 8;
+    end
+end
+
+dom = f.domain(:,1);
+scl = 1./2.^glevel;
+offx = scl*(dom(2)-dom(1));
+offy = scl*(dom(4)-dom(3));
+offz = scl*(dom(6)-dom(5));
+xo = offx.*double(gcol);
+yo = offy.*double(grow);
+zo = offz.*double(gdep);
+gdomain = [dom(1) + xo; dom(1)+offx + xo; ...
+           dom(3) + yo; dom(3)+offy + yo; ...
+           dom(5) + zo; dom(5)+offz + zo];
+
+% Finally, fill in the coefficients
+gid = 1;
+for fid = 1:fnboxes
+    if ( fheight(fid) == 0 )
+        fm = fmorton(fid);
+        %fc = fcol(fid);
+        %fr = frow(fid);
+        fl = flevel(fid);
+        while ( ~(gmorton(gid)==fm && glevel(gid)==fl) && gid <= gnboxes )
+        %while ( ~(gcol(gid)==fc && grow(gid)==fr && glevel(gid)==fl) && gid <= gnboxes )
+            gid = gid + 1;
+        end
+        if ( gid > gnboxes )
+            error('Couldn''t find corresponding tree node.');
+        end
+        gcoeffs{gid} = fcoeffs{fid};
+    end
+end
+
+% Propagate the coefficients to the leaves, since they may be stored in
+% intermediate nodes of g (which used to be leaves of f)
+for l = 1:gnlevels-1
+    for gid = glevelIdx(l):glevelIdx(l+1)-1
+        if ( gheight(gid) > 0 && ~isempty(gcoeffs{gid}) )
+            [LLD, LRD, ULD, URD, LLT, LRT, ULT, URT] = coeffs2children(gcoeffs{gid});
+            gcoeffs{gchildren(1,gid)} = LLD;
+            gcoeffs{gchildren(2,gid)} = LRD;
+            gcoeffs{gchildren(3,gid)} = ULD;
+            gcoeffs{gchildren(4,gid)} = URD;
+            gcoeffs{gchildren(5,gid)} = LLT;
+            gcoeffs{gchildren(6,gid)} = LRT;
+            gcoeffs{gchildren(7,gid)} = ULT;
+            gcoeffs{gchildren(8,gid)} = URT;
+        end
+    end
+end
+
+% Make a new treefun2
+g = treefun3();
+g.n        = f.n;
+g.domain   = gdomain;
+g.level    = glevel;
+g.height   = gheight;
+g.id       = 1:gnboxes;
+g.parent   = gparent;
+g.children = gchildren;
+g.coeffs   = gcoeffs;
+g.col      = gcol;
+g.row      = grow;
+g.morton   = gmorton;
+
+end
+
+%% Column/row routines
+function nbrxy = xycolleagues(xy, level)
+    x = xy(:,1);
+    y = xy(:,2);
+    z = xy(:,3);
+    %bnd = 2^(level-1)-1;
+    bnd = bitshift(uint64(1), level-1) - 1;
+    % Order is important here:
+    %nbrxy = max(uint64(0), min(bnd, [x-1 y ; x+1 y ; x y-1 ; x y+1 ; x-1 y-1 ; x+1 y-1 ; x-1 y+1 ; x+1 y+1]));
+    nbrxy = max(uint64(0), min(bnd, [[x-1 y-1 z-1; x y-1 z-1; x+1 y-1 z-1; x-1 y z-1; x y z-1; x+1 y z-1; x-1 y+1 z-1; x y+1 z-1; x+1 y+1 z-1]; ...
+                                     [x-1 y-1 z  ; x y-1 z  ; x+1 y-1 z  ; x-1 y z  ;          x+1 y z  ; x-1 y+1 z  ; x y+1 z  ; x+1 y+1 z  ]; ...
+                                     [x-1 y-1 z+1; x y-1 z+1; x+1 y-1 z+1; x-1 y z+1; x y z+1; x+1 y z+1; x-1 y+1 z+1; x y+1 z+1; x+1 y+1 z+1]]));
+    %nbrxy = max(uint64(0), min(bnd, [x-1 y-1; x y-1; x+1 y-1 ; x-1 y ; x+1 y ; x-1 y+1 ; x y+1 ; x+1 y+1]));
+end
+
+function pxy = xyparent(xy)
+    pxy = bitshift(xy, -1);
+end
+
+function [cxy1, cxy2, cxy3, cxy4, cxy5, cxy6, cxy7, cxy8] = xychildren(pxy)
+    cxy1 = bitshift(pxy, 2);
+    cxy2 = cxy1 + uint64([1 0 0]);
+    cxy3 = cxy1 + uint64([0 1 0]);
+    cxy4 = cxy1 + uint64([1 1 0]);
+    cxy5 = cxy1 + uint64([0 0 1]);
+    cxy6 = cxy1 + uint64([1 0 1]);
+    cxy7 = cxy1 + uint64([0 1 1]);
+    cxy8 = cxy1 + uint64([1 1 1]);
+end
+
+
+%% Morton routines
+function [xy, idx] = mortonsort(xy)
+    m = xy2morton(xy);
+    [m, idx] = sort(m);
+    xy = morton2xy(m);
+end
+
+function morton = xy2morton(xy)
+    partxy = Part1By2_64(xy);
+    morton = bitshift(partxy(:,3), 2) + bitshift(partxy(:,2), 1) + partxy(:,1);
+end
+
+function xy = morton2xy(morton)
+    xy = Compact1By2_64([morton bitshift(morton, -1) bitshift(morton, -2)]);
+end
+
+function coll = mortoncolleagues(m, level)
+    %bnd = 2^(level-1)-1;
+    bnd = bitshift(uint64(1), level-1) - 1;
+    %[x, y] = morton2cartesian(m);
+    %nbrx = max(uint64(0), min(bnd, [x-1 ; x+1 ; x   ; x   ; x-1 ; x+1 ; x-1 ; x+1]));
+    %nbry = max(uint64(0), min(bnd, [y   ; y   ; y-1 ; y+1 ; y-1 ; y-1 ; y+1 ; y+1]));
+    %coll = cartesian2morton(nbrx, nbry);
+    xy = morton2xy(m);
+    x = xy(:,1);
+    y = xy(:,2);
+    z = xy(:,3);
+    nbrxy = min(bnd, [[x-1 y-1 z-1; x y-1 z-1; x+1 y-1 z-1; x-1 y z-1; x y z-1; x+1 y z-1; x-1 y+1 z-1; x y+1 z-1; x+1 y+1 z-1]; ...
+                      [x-1 y-1 z  ; x y-1 z  ; x+1 y-1 z  ; x-1 y z  ;          x+1 y z  ; x-1 y+1 z  ; x y+1 z  ; x+1 y+1 z  ]; ...
+                      [x-1 y-1 z+1; x y-1 z+1; x+1 y-1 z+1; x-1 y z+1; x y z+1; x+1 y z+1; x-1 y+1 z+1; x y+1 z+1; x+1 y+1 z+1]]);
+    coll = xy2morton(nbrxy);
+end
+
+function p = mortonparent(c)
+    p = bitshift(c, -3);
+end
+
+function [c1, c2, c3, c4, c5, c6, c7, c8] = mortonchildren(p)
+    c1 = bitshift(p, 3);
+    c2 = c1 + 1;
+    c3 = c1 + 2;
+    c4 = c1 + 3;
+    c5 = c1 + 4;
+    c6 = c1 + 5;
+    c7 = c1 + 6;
+    c8 = c1 + 7;
+end
+
+function [x, y, z] = morton2cartesian(morton)
+    x = Compact1By2_64(morton);
+    y = Compact1By2_64(bitshift(morton, -1));
+    z = Compact1By2_64(bitshift(morton, -2));
+end
+
+function morton = cartesian2morton(x, y)
+    morton = bitshift(Part1By2_64(y), 1) + Part1By2_64(x);
+end
+
+function x = Part1By1(x)
+% "Insert" a 0 bit after each of the 16 low bits of x
+    x = bitand(x, 0x0000ffff);
+    x = bitand(bitxor(x, bitshift(x, 8)), 0x00ff00ff);
+    x = bitand(bitxor(x, bitshift(x, 4)), 0x0f0f0f0f);
+    x = bitand(bitxor(x, bitshift(x, 2)), 0x33333333);
+    x = bitand(bitxor(x, bitshift(x, 1)), 0x55555555);
+end
+
+function x = Compact1By1(x)
+% Inverse of Part1By1 - "delete" all odd-indexed bits
+    x = bitand(x, 0x55555555);
+    x = bitand(bitxor(x, bitshift(x, -1)), 0x33333333);
+    x = bitand(bitxor(x, bitshift(x, -2)), 0x0f0f0f0f);
+    x = bitand(bitxor(x, bitshift(x, -4)), 0x00ff00ff);
+    x = bitand(bitxor(x, bitshift(x, -8)), 0x0000ffff);
+end
+
+function x = Part1By2_64(x)
+% "Insert" two 0 bits after each of the 16 low bits of x
+% https://github.com/trevorprater/pymorton/blob/f248683c19d90e904193c58bbd03e77bc2c43768/pymorton/pymorton.py#L20
+
+x = bitand(x, uint64(0x1fffff));
+x = bitand(bitxor(x, bitshift(x, 32)), uint64(0x1f00000000ffff));
+x = bitand(bitxor(x, bitshift(x, 16)), uint64(0x1f0000ff0000ff));
+x = bitand(bitxor(x, bitshift(x, 8)),  uint64(0x100f00f00f00f00f));
+x = bitand(bitxor(x, bitshift(x, 4)),  uint64(0x10c30c30c30c30c3));
+x = bitand(bitxor(x, bitshift(x, 2)),  uint64(0x1249249249249249));
+
+end
+
+function x = Compact1By2_64(x)
+% Inverse of Part1By2 - "delete" all odd-indexed bits
+    x = bitand(x, uint64(0x1249249249249249));
+    x = bitand(bitxor(x, bitshift(x, -2)),  uint64(0x10c30c30c30c30c3));
+    x = bitand(bitxor(x, bitshift(x, -4)),  uint64(0x100f00f00f00f00f));
+    x = bitand(bitxor(x, bitshift(x, -8)),  uint64(0x1f0000ff0000ff));
+    x = bitand(bitxor(x, bitshift(x, -16)), uint64(0x1f00000000ffff));
+    x = bitand(bitxor(x, bitshift(x, -32)), uint64(0x1fffff));
+end
\ No newline at end of file
diff --git a/@treefun3/balancef.m b/@treefun3/balancef.m
new file mode 100644
index 0000000..9c86b8c
--- /dev/null
+++ b/@treefun3/balancef.m
@@ -0,0 +1,824 @@
+function g = balancef(f)
+% fortran vol_tree_fix_lr
+% assume square boxes? therefore scale is computed from dom(1:2)
+% nd needs to be fixed later...
+% flatNeighbors and leafNeighbors also need to be updated in the future
+%
+
+nd = 1; % 
+dpars = zeros(1000,1);
+ipars = zeros(100,1);
+zpars = zeros(10,1);
+
+if ( isempty(f) )
+    g = f;
+    return
+end
+
+% dada
+dom = f.domain(:,1)';
+scale = (dom(2)-dom(1))/2;
+norder = f.n;
+ndim = numel(dom(1:end/2));
+
+% convert to vol_tree_build output format
+iper = 0;
+npbox = norder^ndim;
+grid = zeros(ndim,norder^2);
+nbmax = 2^ndim*f.id(end); % is 2^ndim enough, since some box needs to be refined more than once...
+nlmax = f.height(1); % is this ok?
+centers0 = 1/2*(f.domain(1:2:end,:)+f.domain(2:2:end,:));
+centers = zeros(ndim,nbmax); centers(:,1:size(centers0,2)) = centers0;
+nlevels = f.height(1);
+nboxes = f.id(end);
+boxsize = 2*scale./2.^(0:nlmax); 
+laddr = zeros(2,nlmax+1); 
+for ilev=0:nlmax
+  laddr(1,ilev+1) = sum(f.level<ilev)+1;
+  laddr(2,ilev+1) = sum(f.level<=ilev);
+end
+ilevel = zeros(nbmax,1); ilevel(1:nboxes) = f.level;
+iparent = zeros(nbmax,1); iparent(1:nboxes) = f.parent;
+nchild = zeros(nbmax,1); nchild(1:nboxes) = sum(f.children>0);
+ichild = zeros(2^ndim,nbmax); ichild(:,1:nboxes) = f.children; ichild(ichild==0) = -1; % -1 convention in fortran tree
+
+% this does not belong to fortran tree, but needed for treefun plot
+coeffs = cell(1,nbmax); coeffs(1:nboxes) = f.coeffs(1:nboxes);
+
+% call computecoll to get vol_tree_fix_lr coll info
+nnbors0 = zeros(nboxes,1); nbors0 = zeros(3^ndim,nboxes);
+[nnbors0,nbors0] = ...
+  computecoll(ndim,nlevels,nboxes,laddr,boxsize,...
+                  centers(:,1:nboxes),iparent(1:nboxes),nchild(1:nboxes),ichild(:,1:nboxes),iper,...
+                  nnbors0,nbors0);
+nnbors = zeros(nbmax,1); nbors = zeros(3^ndim,nbmax);
+nnbors(1:nboxes) = nnbors0; nbors(:,1:nboxes) = nbors0;
+
+% now fix lr (with coeffs, which is different from)
+[centers,nboxes,laddr,ilevel,iparent,nchild,ichild,nnbors,nbors,coeffs] = ...
+      vol_tree_fix_lr_treefun(ndim,iper,nd,dpars,zpars,ipars,...
+                               norder,npbox,grid,centers,nlevels,nboxes,boxsize,...
+                               nbmax,nlmax,laddr,ilevel,iparent,nchild,ichild,nnbors,nbors,...
+                               coeffs);
+
+centers = centers(:,1:nboxes);
+ilevel = ilevel(1:nboxes);
+iparent = iparent(1:nboxes);
+nchild = nchild(1:nboxes);
+ichild = ichild(:,1:nboxes);
+nnbors = nnbors(1:nboxes);
+nbors = nbors(:,1:nboxes);
+coeffs = coeffs(1:nboxes);
+
+% convert back to treefun format, try to follow treefun naming convention
+g = f;
+scl                 = boxsize(ilevel+1); % x,y,z all the same
+g.domain            = zeros(2*ndim,nboxes);
+g.domain(1:2:end,:) = centers - 1/2*scl;
+g.domain(2:2:end,:) = centers + 1/2*scl;
+g.level             = ilevel(:)';
+g.id                = 1:nboxes;
+g.parent            = iparent(:)';
+g.children          = ichild;
+g.height            = zeros(1,nboxes);
+g.height(nchild>0)  = 1;
+for k = length(g.id):-1:1
+  if ( ~(g.height(k)==0) )
+    g.height(k) = 1 + max(g.height(g.children(:,k)));
+  end
+end
+% redo some stuff
+g.col               = uint64(zeros(1,nboxes));
+g.row               = uint64(zeros(1,nboxes));
+g.dep               = uint64(zeros(1,nboxes));
+% g.morton            = uint64(zeros(1,nboxes));
+mc = 2^ndim;
+for ilev = 0:nlevels-1
+  for ibox = laddr(1,ilev+1):laddr(2,ilev+1)
+    if nchild(ibox) > 0
+      if ndim == 3 % column row stuff
+        j              = 1;
+        jbox           = g.children(j,ibox);
+        g.col(jbox)    = 2*g.col(ibox);
+        g.row(jbox)    = 2*g.row(ibox);
+        g.dep(jbox)    = 2*g.dep(ibox);
+        % g.morton(jbox) = cartesian2morton(g.col(jbox), g.row(jbox));
+        j              = 2;
+        jbox           = g.children(j,ibox);
+        g.col(jbox)    = 2*g.col(ibox) + 1;
+        g.row(jbox)    = 2*g.row(ibox);
+        g.dep(jbox)    = 2*g.dep(ibox);
+        % g.morton(jbox) = cartesian2morton(g.col(jbox), g.row(jbox));
+        j              = 3;
+        jbox           = g.children(j,ibox);
+        g.col(jbox)    = 2*g.col(ibox);
+        g.row(jbox)    = 2*g.row(ibox) + 1;
+        g.dep(jbox)    = 2*g.dep(ibox);
+        % g.morton(jbox) = cartesian2morton(g.col(jbox), g.row(jbox));
+        j              = 4;
+        jbox           = g.children(j,ibox);
+        g.col(jbox)    = 2*g.col(ibox) + 1;
+        g.row(jbox)    = 2*g.row(ibox) + 1;
+        g.dep(jbox)    = 2*g.dep(ibox);
+        % g.morton(jbox) = cartesian2morton(g.col(jbox), g.row(jbox));
+        j              = 5;
+        jbox           = g.children(j,ibox);
+        g.col(jbox)    = 2*g.col(ibox);
+        g.row(jbox)    = 2*g.row(ibox);
+        g.dep(jbox)    = 2*g.dep(ibox) + 1;
+        
+        j              = 6;
+        jbox           = g.children(j,ibox);
+        g.col(jbox)    = 2*g.col(ibox) + 1;
+        g.row(jbox)    = 2*g.row(ibox);
+        g.dep(jbox)    = 2*g.dep(ibox) + 1;
+        
+        j              = 7;
+        jbox           = g.children(j,ibox);
+        g.col(jbox)    = 2*g.col(ibox);
+        g.row(jbox)    = 2*g.row(ibox) + 1;
+        g.dep(jbox)    = 2*g.dep(ibox) + 1;
+        
+        j              = 8;
+        jbox           = g.children(j,ibox);
+        g.col(jbox)    = 2*g.col(ibox) + 1;
+        g.row(jbox)    = 2*g.row(ibox) + 1;
+        g.dep(jbox)    = 2*g.dep(ibox) + 1;
+        
+      end
+    end
+  end
+end
+g.coeffs = coeffs;
+end
+
+function [centers,nboxes,laddr,ilevel,iparent,nchild,ichild,nnbors,nbors,coeffs] = vol_tree_fix_lr_treefun(ndim,iperiod,nd,dpars,zpars,ipars, ...
+                                                                                    norder,npbox,grid,centers,nlevels,nboxes,boxsize, ...
+                                                                                    nbmax,nlmax,laddr,ilevel,iparent,nchild,ichild,nnbors,nbors, ...
+                                                                                    coeffs)
+
+laddrtail=zeros(2,nlmax+1);
+% boxsize(0:nlmax),laddr(2,0:nlmax),laddrtail(2,0:nlmax)
+idx1 = 1;
+%
+bs0=boxsize(0+idx1);
+%
+mc=2^ndim;
+mnbors=3^ndim;
+
+iflag = zeros(nbmax,1); iflag2 = zeros(nbmax,1);
+
+for i=1:nboxes
+   iflag(i) = 0;
+end
+
+% Flag boxes that violate level restriction by "1" Violatioin refers to any
+% box that is directly touching a box that is more than one level finer
+% Method:
+% 1) Carry out upward pass. For each box B, look at the colleagues of B's
+% grandparent 
+% 2) See if any of those colleagues are childless and in contact with B.
+%  
+% Note that we only need to get up to level two, as we will not find a
+% violation at level 0 and level 1 For such boxes, we set iflag(i) = 1
+% figure(1), clf,
+for ilev=nlevels:-1:2
+  distest = 1.05d0*(boxsize(ilev-1+idx1) + boxsize(ilev-2+idx1))/2.0d0;
+  for ibox = laddr(1,ilev+idx1):laddr(2,ilev+idx1) % look at current box
+    idad = iparent(ibox); % its parent box
+    igranddad = iparent(idad); % its grand parent box
+    
+    for i=1:nnbors(igranddad) % look at neighbor of grand parent
+       jbox = nbors(i,igranddad);
+       if((nchild(jbox)==0) && (iflag(jbox)==0)) % if neighbor of grand parent has no child, and has not been flagged
+          ict = 0; % count overlapping dimension between ibox & its grand parent's neighbor
+          for k=1:ndim
+             dis = centers(k,jbox) - centers(k,idad);
+             if (iperiod==1)
+                dp1=dis+bs0;
+                dm1=dis-bs0;
+                if (abs(dis).gt.abs(dp1)), dis=dp1; end
+                if (abs(dis).gt.abs(dm1)), dis=dm1; end
+             end
+             if(abs(dis)<=distest), ict = ict + 1; end
+          end
+          if(ict==ndim)
+            iflag(jbox) = 1;
+            % figure(1),
+            % plotibox(centers(:,ibox),boxsize(ilev-1+idx1)/2);
+            % plotibox(centers(:,jbox),4*boxsize(ilev-1+idx1)/2);
+            % axis equal
+          end
+       end
+    end
+  end
+end
+
+% Find all boxes that need to be given a flag+ A flag+ box will be denoted
+% by setting iflag(box) = 2 This refers to any box that is not already
+% flagged and is bigger than and is contacting a flagged box or another box
+% that has already been given a flag +. It is found by performing an upward
+% pass and looking at the flagged box's parents colleagues and a flag+
+% box's parents colleagues and seeing if they are childless and present the
+% case where a bigger box is contacting a flagged or flag+ box.
+% figure(2), clf,
+%% if secondary violator, 1st send it downstairs during Step 4, then wait to be processed in Step 5 or no need to process at all, seems to be later, from the "worst-case" example 
+for ilev = nlevels:-1:1
+  distest = 1.05d0*(boxsize(ilev+idx1) + boxsize(ilev-1+idx1))/2.0d0;
+  for ibox = laddr(1,ilev+idx1):laddr(2,ilev+idx1)
+    if((iflag(ibox)==1) || (iflag(ibox)==2)) % look at current box
+      idad = iparent(ibox); % its parent box
+      for i=1:nnbors(idad) % look at neighbor of parent
+        jbox = nbors(i,idad);
+        if((nchild(jbox)==0) && (iflag(jbox)==0)) % if neighbor of parent has no child, and has not been flagged
+          ict = 0; % count overlapping dimension between ibox & its parent's neighbor
+          for k=1:ndim
+            dis = centers(k,jbox) - centers(k,ibox);
+            if (iperiod==1)
+              dp1=dis+bs0;
+              dm1=dis-bs0;
+              if (abs(dis)>abs(dp1)), dis=dp1; end
+              if (abs(dis)>abs(dm1)), dis=dm1; end
+            end
+            if(abs(dis)<=distest), ict = ict + 1; end
+          end
+          if(ict==ndim)
+            iflag(jbox) = 2;
+            % figure(2),
+            % plotibox(centers(:,ibox),boxsize(ilev-1+idx1)/2);
+            % plotibox(centers(:,jbox),4*boxsize(ilev-1+idx1)/2);
+            % axis equal
+          end
+        end
+      end
+    end
+  end
+end
+
+% Subdivide all flag and flag+ boxes. Flag all the children of flagged
+% boxes as flag++. Flag++ boxes are denoted by setting iflag(box) = 3. The
+% flag++ boxes need to be checked later to see which of them need further
+% refinement. While creating new boxes, we will need to update all the tree
+% structures as well. Note that all the flagged boxes live between levels 1
+% and nlevels - 2. We process the boxes via a downward pass. We first
+% determine the number of boxes that are going to be subdivided at each
+% level and everything else accordingly
+for ilev = 0:nlevels
+   laddrtail(1,ilev+idx1) = 0;
+   laddrtail(2,ilev+idx1) = -1;
+end
+
+nboxes0 = nboxes;
+%% Frank's thesis, Appendix A, Step 4, divide all of the primary and secondary violators once. I think
+% figure out iflag = 3 case
+for ilev = 1:nlevels-2
+
+  laddrtail(1,ilev+1+idx1) = nboxes+1;
+
+  nbloc = laddr(2,ilev+idx1)-laddr(1,ilev+idx1)+1;
+  
+  [iflag,centers,nboxes,ilevel,iparent,nchild,ichild,coeffs]  = vol_tree_refine_boxes_flag_treefun(ndim,iflag,nd,npbox,...
+                                   dpars,zpars,ipars,grid,nbmax,laddr(1,ilev+idx1),nbloc, ...
+                                   centers,boxsize(ilev+1+idx1),nboxes,ilev,ilevel,iparent,nchild,ichild, ...
+                                   coeffs);
+  
+  laddrtail(2,ilev+1+idx1) = nboxes;
+
+end
+
+ishuffle1 = zeros(nboxes,1);
+centerstmp = centers(:,1:nboxes); ileveltmp = ilevel(1:nboxes); iparenttmp = iparent(1:nboxes); nchildtmp = nchild(1:nboxes); ichildtmp = ichild(:,1:nboxes); iflagtmp = iflag(1:nboxes);
+coeffstmp = coeffs(1:nboxes);
+[centerstmp,laddr,ileveltmp,ishuffle1,iparenttmp,nchildtmp,ichildtmp,coeffstmp,iflagtmp] = ... 
+              vol_tree_reorg_treefun(ndim,nboxes,nd,npbox,centerstmp,nlevels,laddr,...
+                                      laddrtail,ileveltmp,ishuffle1,iparenttmp,nchildtmp,ichildtmp,coeffstmp,iflagtmp);
+centers(:,1:nboxes) = centerstmp; ilevel(1:nboxes) = ileveltmp; iparent(1:nboxes) = iparenttmp; nchild(1:nboxes) = nchildtmp; ichild(:,1:nboxes) = ichildtmp; iflag(1:nboxes) = iflagtmp;
+coeffs(1:nboxes) = coeffstmp;
+
+nnborstmp = nnbors(1:nboxes); nborstmp = nbors(:,1:nboxes);
+[nnborstmp,nborstmp] = computecoll(ndim,nlevels,nboxes,laddr,boxsize,...
+                                      centers(:,1:nboxes),iparent(1:nboxes),nchild(1:nboxes),...
+                                      ichild(:,1:nboxes),iperiod,nnborstmp,nborstmp);
+nnbors(1:nboxes) = nnborstmp; nbors(:,1:nboxes) = nborstmp;
+
+
+% Processing of flag and flag+ boxes is done Start processing flag++ boxes.
+% We will use a similar strategy as before. We keep checking the flag++
+% boxes that require subdivision if they still violate the level
+% restriction criterion, create the new boxes, append them to the end of
+% the list to begin with and in the end reorganize the tree structure. We
+% shall accomplish this via a downward pass as new boxes that get added in
+% the downward pass will also be processed simultaneously. We shall
+% additionally also need to keep on updating the colleague information as
+% we proceed in the downward pass
+% 
+% Reset the flags array to remove all the flag and flag+ cases. This is to
+% ensure reusability of the subdivide_flag routine to handle the flag++ case
+
+for ibox=1:nboxes
+  if(iflag(ibox)~=3), iflag(ibox) = 0; end
+end
+
+for ilev = 0:nlevels
+  laddrtail(1,ilev+idx1) = 0;
+  laddrtail(2,ilev+idx1) = -1;
+end
+
+% boxsize(0:nlmax),laddr(2,0:nlmax),laddrtail(2,0:nlmax)
+%% Frank's thesis, Appendix A, Step 5, at each level, for a descendant of a primary violator, test whether it is still in violation
+for ilev = 2:nlevels-2
+
+  % Step 1
+  iflagtmp = iflag(1:nboxes);
+  iflagtmp = vol_updateflags(ndim,iperiod,ilev,nboxes,nlevels, ...
+            laddr,nchild(1:nboxes),ichild(:,1:nboxes),nnbors(1:nboxes),nbors(:,1:nboxes),centers(:,1:nboxes),boxsize,iflagtmp);
+  iflag(1:nboxes) = iflagtmp;
+
+  iflagtmp = iflag(1:nboxes);
+  iflagtmp = vol_updateflags(ndim,iperiod,ilev,nboxes,nlevels, ...
+            laddrtail,nchild(1:nboxes),ichild(:,1:nboxes),nnbors(1:nboxes),nbors(:,1:nboxes),centers(:,1:nboxes),boxsize,iflagtmp);
+  iflag(1:nboxes) = iflagtmp;
+
+  % Step 2
+  laddrtail(1,ilev+1+idx1) = nboxes + 1;
+
+  nbloc = laddr(2,ilev+idx1)-laddr(1,ilev+idx1)+1;
+  [iflag,centers,nboxes,ilevel,iparent,nchild,ichild,coeffs] = vol_tree_refine_boxes_flag_treefun(ndim,iflag,nd,npbox, ...
+                                                       dpars,zpars,ipars,grid,nbmax,laddr(1,ilev+idx1),nbloc, ...
+                                                       centers,boxsize(ilev+1+idx1),nboxes,ilev,ilevel,iparent,nchild,ichild, ...
+                                                       coeffs);
+  
+  nbloc = laddrtail(2,ilev+idx1)-laddrtail(1,ilev+idx1)+1;
+  [iflag,centers,nboxes,ilevel,iparent,nchild,ichild,coeffs] = vol_tree_refine_boxes_flag_treefun(ndim,iflag,nd,npbox, ...
+                                                       dpars,zpars,ipars,grid,nbmax,laddrtail(1,ilev+idx1),nbloc, ...
+                                                       centers,boxsize(ilev+1+idx1),nboxes,ilev,ilevel,iparent,nchild,ichild, ...
+                                                       coeffs);
+  laddrtail(2,ilev+1+idx1) = nboxes;      
+
+  % Step 3
+  for ibox = laddrtail(1,ilev+1+idx1):laddrtail(2,ilev+1+idx1)
+    nnbors(ibox) = 0;
+    idad = iparent(ibox);
+    for i=1:nnbors(idad)
+      jbox = nbors(i,idad);
+      for j=1:mc
+        kbox = ichild(j,jbox);
+        if(kbox>0) 
+          ifnbor=1;
+          for k=1:ndim
+            dis=centers(k,kbox)-centers(k,ibox);
+            if (iperiod==1) 
+              dp1=dis+bs0;
+              dm1=dis-bs0;
+              if (abs(dis)>abs(dp1)), dis=dp1; end
+              if (abs(dis)>abs(dm1)), dis=dm1; end
+            end
+            if((abs(dis)>1.05*boxsize(ilev+1+idx1)))
+              ifnbor=0;
+              break
+            end
+          end
+          if (ifnbor==1)
+            nnbors(ibox) = nnbors(ibox)+1;
+            nbors(nnbors(ibox),ibox) = kbox;
+          end
+        end
+      end
+    end
+    % End of computing colleagues of box i
+  end
+  
+end
+
+ishuffle2 = zeros(nboxes,1);
+centerstmp = centers(:,1:nboxes); ileveltmp = ilevel(1:nboxes); iparenttmp = iparent(1:nboxes); nchildtmp = nchild(1:nboxes); ichildtmp = ichild(:,1:nboxes); iflagtmp = iflag(1:nboxes);
+coeffstmp = coeffs(1:nboxes);
+[centerstmp,laddr,ileveltmp,ishuffle2,iparenttmp,nchildtmp,ichildtmp,coeffstmp,iflagtmp] = ... 
+              vol_tree_reorg_treefun(ndim,nboxes,nd,npbox,centerstmp,nlevels,laddr,...
+                                      laddrtail,ileveltmp,ishuffle2,iparenttmp,nchildtmp,ichildtmp,coeffstmp,iflagtmp);
+centers(:,1:nboxes) = centerstmp; ilevel(1:nboxes) = ileveltmp; iparent(1:nboxes) = iparenttmp; nchild(1:nboxes) = nchildtmp; ichild(:,1:nboxes) = ichildtmp; iflag(1:nboxes) = iflagtmp;
+coeffs(1:nboxes) = coeffstmp;
+
+nnborstmp = nnbors(1:nboxes); nborstmp = nbors(:,1:nboxes);
+[nnborstmp,nborstmp] = computecoll(ndim,nlevels,nboxes,laddr,boxsize,...
+                                      centers(:,1:nboxes),iparent(1:nboxes),nchild(1:nboxes),...
+                                      ichild(:,1:nboxes),iperiod,nnborstmp,nborstmp);
+nnbors(1:nboxes) = nnborstmp; nbors(:,1:nboxes) = nborstmp;
+
+end
+
+function [iflag,centers,nbctr,ilevel,iparent,nchild,ichild,coeffs] = vol_tree_refine_boxes_flag_treefun(ndim,iflag,nd,npbox,dpars,zpars,ipars,grid,nboxes,ifirstbox,nbloc,centers,bs,nbctr,nlctr,ilevel,iparent,nchild,ichild,coeffs)
+
+% myfun = @(x,y,z) sin(x).*cos(y).*sin(z)+exp(x.*y.*z);
+% myvals = myfun(2*(xx0-1/2),2*(yy0-1/2),2*(zz0-1/2));
+% mycoeffs = treefun3.vals2coeffs(myvals);
+% mycvals1_2 = myfun(cxx1,cyy1,czz1);
+% mycvals1 = reshape(treefun3.coeffs2checkvals(mycoeffs,cxx1(:),cyy1(:),czz1(:)),[norder norder norder]);
+% mycvals2 = reshape(treefun3.coeffs2checkvals(mycoeffs,cxx2(:),cyy2(:),czz2(:)),[norder norder norder]);
+% mycvals3 = reshape(treefun3.coeffs2checkvals(mycoeffs,cxx3(:),cyy3(:),czz3(:)),[norder norder norder]);
+% mycvals4 = reshape(treefun3.coeffs2checkvals(mycoeffs,cxx4(:),cyy4(:),cyy4(:)),[norder norder norder]);
+% mycvals5 = reshape(treefun3.coeffs2checkvals(mycoeffs,cxx5(:),cyy5(:),czz5(:)),[norder norder norder]);
+% mycvals6 = reshape(treefun3.coeffs2checkvals(mycoeffs,cxx6(:),cyy6(:),czz6(:)),[norder norder norder]);
+% mycvals7 = reshape(treefun3.coeffs2checkvals(mycoeffs,cxx7(:),cyy7(:),czz7(:)),[norder norder norder]);
+% mycvals8 = reshape(treefun3.coeffs2checkvals(mycoeffs,cxx8(:),cyy8(:),cyy8(:)),[norder norder norder]);
+% 
+% figure(2), clf, scatter3(2*(xx0(:)-1/2),2*(yy0(:)-1/2),2*(zz0(:)-1/2),[],log10(abs(vals(:))))
+% hold on, scatter3(cxx1(:),cyy1(:),czz1(:),[],log10(abs(cvals1(:))))
+% hold on, scatter3(cxx2(:),cyy2(:),czz2(:),[],log10(abs(cvals2(:))))
+% hold on, scatter3(cxx3(:),cyy3(:),czz3(:),[],log10(abs(cvals3(:))))
+% hold on, scatter3(cxx4(:),cyy4(:),czz4(:),[],log10(abs(cvals4(:))))
+% hold on, scatter3(cxx5(:),cyy5(:),czz5(:),[],log10(abs(cvals5(:))))
+% hold on, scatter3(cxx6(:),cyy6(:),czz6(:),[],log10(abs(cvals6(:))))
+% hold on, scatter3(cxx7(:),cyy7(:),czz7(:),[],log10(abs(cvals7(:))))
+% hold on, scatter3(cxx8(:),cyy8(:),czz8(:),[],log10(abs(cvals8(:))))
+
+% additional
+persistent xx0 yy0 xxx0 yyy0 zz0 zzz0 nstored
+norder=round((npbox)^(1/ndim));
+n=norder;
+
+if ndim == 3
+  nalias = n;
+  nrefpts = n; % Sample at equispaced points to test error
+
+  if ( isempty(xx0) || isempty(xxx0) || n ~= nstored )
+      nstored = n;
+      x0 = (1-cos(pi*(2*(1:nalias)'-1)/(2*nalias)))/2;
+      [xx0, yy0, zz0] = ndgrid(x0);
+      [xxx0, yyy0, zzz0] = ndgrid(linspace(0, 1, nrefpts));
+  end
+  sclx = 2*bs;
+  scly = 2*bs;
+  sclz = 2*bs;
+end
+
+%
+isgn=zeros(ndim,2^ndim);
+isgn=get_child_box_sign(ndim,isgn);
+
+mc=2^ndim;
+
+ilastbox = ifirstbox+nbloc-1;
+
+bsh = bs/2;
+
+isum=zeros(nbloc,1);
+
+isum = cumsum(iflag(ifirstbox:(ifirstbox+nbloc-1))>0); % there are iflag = 2 entries
+
+for i = 1:nbloc
+  ibox = ifirstbox + i-1;
+  if(iflag(ibox)>0)
+    nbl = nbctr + (isum(i)-1)*mc;
+    nchild(ibox) = mc;
+    % additional
+    if ndim == 3
+      xx = sclx*xx0 + (centers(1,ibox)-1/2*sclx);
+      yy = scly*yy0 + (centers(2,ibox)-1/2*scly);
+      zz = sclz*zz0 + (centers(3,ibox)-1/2*sclz);
+      vals = treefun3.coeffs2vals(coeffs{ibox});
+      % split into eight child boxes
+      dom = [-1 1 -1 1 -1 1];
+      xmid = mean(dom(1:2));
+      ymid = mean(dom(3:4));
+      zmid = mean(dom(5:6));
+      j=1;
+      jbox = nbl+j;
+      cdom1              = [dom(1) xmid dom(3) ymid dom(5) zmid];
+      csclx1             = diff(cdom1(1:2));
+      cscly1             = diff(cdom1(3:4));
+      csclz1             = diff(cdom1(5:6));
+      cxx1               = csclx1*xx0 + cdom1(1); 
+      cyy1               = cscly1*yy0 + cdom1(3); 
+      czz1               = csclz1*zz0 + cdom1(5); 
+      cvals1             = reshape(treefun3.coeffs2checkvals(coeffs{ibox},cxx1(:),cyy1(:),czz1(:)),[norder norder norder]);
+      ccoeffs1           = treefun3.vals2coeffs(cvals1);
+      coeffs{jbox}       = ccoeffs1;
+      j=2;
+      jbox = nbl+j;
+      cdom2              = [xmid dom(2) dom(3) ymid dom(5) zmid];
+      csclx2             = diff(cdom2(1:2));
+      cscly2             = diff(cdom2(3:4));
+      csclz2             = diff(cdom2(5:6));
+      cxx2               = csclx2*xx0 + cdom2(1); 
+      cyy2               = cscly2*yy0 + cdom2(3); 
+      czz2               = csclz2*zz0 + cdom2(5); 
+      cvals2             = reshape(treefun3.coeffs2checkvals(coeffs{ibox},cxx2(:),cyy2(:),czz2(:)),[norder norder norder]);
+      ccoeffs2           = treefun3.vals2coeffs(cvals2);
+      coeffs{jbox}       = ccoeffs2;
+      j=3;
+      jbox = nbl+j;
+      cdom3              = [dom(1) xmid ymid dom(4) dom(5) zmid];
+      csclx3             = diff(cdom3(1:2));
+      cscly3             = diff(cdom3(3:4));
+      csclz3             = diff(cdom3(5:6));
+      cxx3               = csclx3*xx0 + cdom3(1); 
+      cyy3               = cscly3*yy0 + cdom3(3); 
+      czz3               = csclz3*zz0 + cdom3(5); 
+      cvals3             = reshape(treefun3.coeffs2checkvals(coeffs{ibox},cxx3(:),cyy3(:),czz3(:)),[norder norder norder]);
+      ccoeffs3           = treefun3.vals2coeffs(cvals3);
+      coeffs{jbox}       = ccoeffs3;
+      j=4;
+      jbox = nbl+j;
+      cdom4              = [xmid dom(2) ymid dom(4) dom(5) zmid];
+      csclx4             = diff(cdom4(1:2));
+      cscly4             = diff(cdom4(3:4));
+      csclz4             = diff(cdom4(5:6));
+      cxx4               = csclx4*xx0 + cdom4(1); 
+      cyy4               = cscly4*yy0 + cdom4(3); 
+      czz4               = csclz4*zz0 + cdom4(5); 
+      cvals4             = reshape(treefun3.coeffs2checkvals(coeffs{ibox},cxx4(:),cyy4(:),czz4(:)),[norder norder norder]);
+      ccoeffs4           = treefun3.vals2coeffs(cvals4);
+      coeffs{jbox}       = ccoeffs4;
+      j=5;
+      jbox = nbl+j;
+      cdom5              = [dom(1) xmid dom(3) ymid zmid dom(6)];
+      csclx5             = diff(cdom5(1:2));
+      cscly5             = diff(cdom5(3:4));
+      csclz5             = diff(cdom5(5:6));
+      cxx5               = csclx5*xx0 + cdom5(1); 
+      cyy5               = cscly5*yy0 + cdom5(3); 
+      czz5               = csclz5*zz0 + cdom5(5); 
+      cvals5             = reshape(treefun3.coeffs2checkvals(coeffs{ibox},cxx5(:),cyy5(:),czz5(:)),[norder norder norder]);
+      ccoeffs5           = treefun3.vals2coeffs(cvals5);
+      coeffs{jbox}       = ccoeffs5;
+      j=6;
+      jbox = nbl+j;
+      cdom6              = [xmid dom(2) dom(3) ymid zmid dom(6)];
+      csclx6             = diff(cdom6(1:2));
+      cscly6             = diff(cdom6(3:4));
+      csclz6             = diff(cdom6(5:6));
+      cxx6               = csclx6*xx0 + cdom6(1); 
+      cyy6               = cscly6*yy0 + cdom6(3); 
+      czz6               = csclz6*zz0 + cdom6(5); 
+      cvals6             = reshape(treefun3.coeffs2checkvals(coeffs{ibox},cxx6(:),cyy6(:),czz6(:)),[norder norder norder]);
+      ccoeffs6           = treefun3.vals2coeffs(cvals6);
+      coeffs{jbox}       = ccoeffs6;
+      j=7;
+      jbox = nbl+j;
+      cdom7              = [dom(1) xmid ymid dom(4) zmid dom(6)];
+      csclx7             = diff(cdom7(1:2));
+      cscly7             = diff(cdom7(3:4));
+      csclz7             = diff(cdom7(5:6));
+      cxx7               = csclx7*xx0 + cdom7(1); 
+      cyy7               = cscly7*yy0 + cdom7(3); 
+      czz7               = csclz7*zz0 + cdom7(5); 
+      cvals7             = reshape(treefun3.coeffs2checkvals(coeffs{ibox},cxx7(:),cyy7(:),czz7(:)),[norder norder norder]);
+      ccoeffs7           = treefun3.vals2coeffs(cvals7);
+      coeffs{jbox}       = ccoeffs7;
+      j=8;
+      jbox = nbl+j;
+      cdom8              = [xmid dom(2) ymid dom(4) zmid dom(6)];
+      csclx8             = diff(cdom8(1:2));
+      cscly8             = diff(cdom8(3:4));
+      csclz8             = diff(cdom8(5:6));
+      cxx8               = csclx8*xx0 + cdom8(1); 
+      cyy8               = cscly8*yy0 + cdom8(3); 
+      czz8               = csclz8*zz0 + cdom8(5); 
+      cvals8             = reshape(treefun3.coeffs2checkvals(coeffs{ibox},cxx8(:),cyy8(:),czz8(:)),[norder norder norder]);
+      ccoeffs8           = treefun3.vals2coeffs(cvals8);
+      coeffs{jbox}       = ccoeffs8;
+      %
+      coeffs{ibox}       = [];
+    end
+    for j=1:mc
+      jbox = nbl+j;
+      for k=1:ndim
+        centers(k,jbox) = centers(k,ibox)+isgn(k,j)*bsh;
+      end
+      iparent(jbox) = ibox;
+      nchild(jbox) = 0;
+      for l=1:mc
+        ichild(l,jbox) = -1;
+      end
+      ichild(j,ibox) = jbox;
+      ilevel(jbox) = nlctr+1; 
+      if(iflag(ibox)==1), iflag(jbox) = 3; end
+      if(iflag(ibox)==2), iflag(jbox) = 0; end
+    end
+  end
+end
+
+if(nbloc>0), nbctr = nbctr + isum(nbloc)*mc; end
+
+end
+
+function iflag = vol_updateflags(ndim,iperiod,curlev,nboxes,nlevels,laddr,nchild,ichild,nnbors,nbors,centers,boxsize,iflag)
+%
+%
+
+idx1 = 1;
+bs0=boxsize(0+idx1);
+
+mc=2^ndim;
+distest = 1.05d0*(boxsize(curlev+idx1) + boxsize(curlev+1+idx1))/2.0d0;
+
+for ibox = laddr(1,curlev+idx1):laddr(2,curlev+idx1)
+  if(iflag(ibox)==3)
+    iflag(ibox) = 0;
+    for i=1:nnbors(ibox)
+      jbox = nbors(i,ibox);
+      
+      for j=1:mc
+        kbox = ichild(j,jbox);
+        if(kbox>0)
+          if(nchild(kbox)>0)
+            ict = 0;
+            for k=1:ndim
+              dis = centers(k,kbox) - centers(k,ibox);
+              if (iperiod==1) 
+                dp1=dis+bs0;
+                dm1=dis-bs0;
+                if (abs(dis)>abs(dp1)), dis=dp1; end
+                if (abs(dis)>abs(dm1)), dis=dm1; end
+              end
+              if(abs(dis)<=distest), ict = ict + 1; end
+            end
+            if(ict==ndim)
+               iflag(ibox) = 1;
+               break;
+            end
+          end
+        end
+      end
+      if iflag(ibox) == 1
+       break;
+      end
+    end
+    % 1111       continue   
+  end
+end
+
+end
+
+function [nnbors,nbors] = computecoll(ndim,nlevels,nboxes,laddr,boxsize,centers,iparent,nchild,ichild,iperiod,nnbors,nbors)
+
+% laddr(2,0:nlevels)
+% boxsize(0:nlevels)
+
+idx1 = 1;
+
+bs0=boxsize(0+idx1);
+
+mc= 2^ndim;
+mnbors=3^ndim;
+
+for i=1:nboxes
+  nnbors(i) = 0;
+  for j=1:mnbors
+    nbors(j,i) = -1;
+  end
+end
+
+nnbors(1) = 1;
+nbors(1,1) = 1;
+for ilev = 1:nlevels
+  ifirstbox = laddr(1,ilev+idx1);
+  ilastbox = laddr(2,ilev+idx1);
+  
+  for ibox = ifirstbox:ilastbox
+    dad = iparent(ibox);
+    for i=1:nnbors(dad)
+      jbox = nbors(i,dad);
+      for j=1:mc
+        kbox = ichild(j,jbox);
+        if(kbox>0)
+          ifnbor=1;
+          for k=1:ndim
+            dis=abs(centers(k,kbox)-centers(k,ibox));
+            if (iperiod==1)
+              dp1=bs0-dis;
+              if (dp1<dis), dis=dp1; end
+            end
+            if (dis>1.05*boxsize(ilev+idx1))
+              ifnbor=0;
+              break
+            end
+          end
+             
+          if(ifnbor==1)
+            nnbors(ibox) = nnbors(ibox)+1;
+            nbors(nnbors(ibox),ibox) = kbox;
+          end
+        end
+      end
+    end
+  end
+end
+
+end
+
+function [centers,laddr,ilevel,ishuffle,iparent,nchild,ichild,coeffs,iflag] = vol_tree_reorg_treefun(ndim,nboxes,nd,npbox,centers,nlevels,laddr,laddrtail,ilevel,ishuffle,iparent,nchild,ichild,coeffs,iflag)
+%
+%
+
+% laddrtail(2,0:nlevels)
+% laddr(2,0:nlevels)
+% tladdr(2,0:nlevels)
+
+tladdr = zeros(2,nlevels+1);
+
+idx1 = 1;
+
+mc=2^ndim;
+
+tilevel=zeros(nboxes,1); tiparent=zeros(nboxes,1); tnchild=zeros(nboxes,1);
+tichild=zeros(mc,nboxes); tiflag=zeros(nboxes,1); iboxtocurbox=zeros(nboxes,1);
+tcenters=zeros(ndim,nboxes);
+
+for ilev = 0:nlevels
+  tladdr(1,ilev+idx1) = laddr(1,ilev+idx1);
+  tladdr(2,ilev+idx1) = laddr(2,ilev+idx1);
+end
+
+tcenters = centers;
+tilevel = ilevel;
+tiparent = iparent;
+tnchild = nchild;
+tichild = ichild;
+tcoeffs = coeffs;
+
+for ibox=1:nboxes
+  tiflag(ibox) = iflag(ibox);
+end
+
+for ilev = 0:1
+  for ibox = laddr(1,ilev+idx1):laddr(2,ilev+idx1)
+    iboxtocurbox(ibox) = ibox;
+  end
+end
+
+ilevptr=zeros(nlevels+1,1); ilevptr2=zeros(nlevels,1);
+
+ilevptr(2) = laddr(1,2+idx1);
+
+for ilev=2:nlevels
+  nblev = laddr(2,ilev+idx1)-laddr(1,ilev+idx1)+1; % original num of boxes at ilev
+  ilevptr2(ilev) = ilevptr(ilev) + nblev;
+  nblev = laddrtail(2,ilev+idx1)-laddrtail(1,ilev+idx1)+1; % additional number of boxes at ilev due to refinement
+  ilevptr(ilev+1) = ilevptr2(ilev) + nblev;
+end
+
+curbox = laddr(1,2+idx1);
+ishuffle = zeros(nboxes,1);
+for ilev=2:nlevels
+  laddr(1,ilev+idx1) = curbox;
+  for ibox = tladdr(1,ilev+idx1):tladdr(2,ilev+idx1)
+    ilevel(curbox) = tilevel(ibox);
+    nchild(curbox) = tnchild(ibox);
+    ishuffle(ibox) = curbox; % ibox info copied to curbox
+    for k=1:ndim
+      centers(k,curbox) = tcenters(k,ibox);
+    end
+    coeffs{curbox} = tcoeffs{ibox};
+    iflag(curbox) = tiflag(ibox);
+    iboxtocurbox(ibox) = curbox;
+    
+    curbox = curbox + 1;
+  end
+  for ibox = laddrtail(1,ilev+idx1):laddrtail(2,ilev+idx1)
+    ilevel(curbox) = tilevel(ibox);
+    ishuffle(ibox) = curbox; % ibox info copied to curbox
+    for k=1:ndim
+      centers(k,curbox) = tcenters(k,ibox);
+    end
+    coeffs{curbox} = tcoeffs{ibox};
+    nchild(curbox) = tnchild(ibox);
+    iflag(curbox) = tiflag(ibox);
+    iboxtocurbox(ibox) = curbox;
+    
+    curbox = curbox + 1;
+  end
+  laddr(2,ilev+idx1) = curbox-1;
+end
+
+for ibox=1:nboxes
+  if(tiparent(ibox)==-1), iparent(iboxtocurbox(ibox)) = -1; end
+  if(tiparent(ibox)>0)
+    iparent(iboxtocurbox(ibox)) = iboxtocurbox(tiparent(ibox)); end
+  for i=1:mc
+    if(tichild(i,ibox)==-1), ichild(i,iboxtocurbox(ibox)) = -1; end
+    if(tichild(i,ibox)>0)
+      ichild(i,iboxtocurbox(ibox)) = iboxtocurbox(tichild(i,ibox)); end
+  end
+end
+
+end
+
+function isgn = get_child_box_sign(ndim,isgn)
+
+mc = 2^ndim;
+for j=1:ndim
+  isgn(j,1)=-1;
+end
+
+for j=1:ndim
+  for i=1:2^(j-1):mc
+    if (i>1), isgn(j,i)=-isgn(j,i-2^(j-1)); end
+    for k=1:2^(j-1)-1
+      isgn(j,i+k)=isgn(j,i);
+    end
+  end
+end
+
+end
\ No newline at end of file
diff --git a/@treefun3/coeffs2checkvals.m b/@treefun3/coeffs2checkvals.m
new file mode 100644
index 0000000..14014a3
--- /dev/null
+++ b/@treefun3/coeffs2checkvals.m
@@ -0,0 +1,28 @@
+function vals = coeffs2checkvals(coeffs,x,y,z)
+%
+%
+
+persistent Evalx Evaly Evalz
+[p,~,~,nd] = size(coeffs);
+ncheckpts = numel(x); % hopefully not too large...
+
+Evalx = ones(ncheckpts, p);
+Evaly = ones(ncheckpts, p);
+Evalz = ones(ncheckpts, p);
+Evalx(:,2) = x(:);
+Evaly(:,2) = y(:);
+Evalz(:,2) = z(:);
+for k=3:p
+  Evalx(:,k) = 2*x(:).*Evalx(:,k-1)-Evalx(:,k-2);
+  Evaly(:,k) = 2*y(:).*Evaly(:,k-1)-Evaly(:,k-2);
+  Evalz(:,k) = 2*z(:).*Evalz(:,k-1)-Evalz(:,k-2);
+end
+% coeffs to value map
+vals = zeros(nd,ncheckpts);
+for k=1:ncheckpts
+  tmp1 = permute(tensorprod(Evalx(k,:),coeffs,2,1),[2 3 1 4]);
+  tmp2 = permute(tensorprod(Evaly(k,:),tmp1,2,1),[2 3 1 4]);
+  vals(:,k) = squeeze(permute(tensorprod(Evalz(k,:),tmp2,2,1),[2 3 1 4]));
+end
+
+end
\ No newline at end of file
diff --git a/@treefun3/coeffs2refvals.m b/@treefun3/coeffs2refvals.m
new file mode 100644
index 0000000..080133b
--- /dev/null
+++ b/@treefun3/coeffs2refvals.m
@@ -0,0 +1,23 @@
+function vals = coeffs2refvals(coeffs)
+%
+%
+
+persistent Eval pstored
+[p,~,~,nd] = size(coeffs);
+nrefpts = p;
+
+if ( isempty(Eval) || p ~= pstored )
+    pstored = p;
+    Eval = ones(nrefpts, p);
+    x = linspace(-1, 1, nrefpts).';
+    Eval(:,2) = x;
+    for k=3:p
+        Eval(:,k) = 2*x.*Eval(:,k-1)-Eval(:,k-2);
+    end
+end
+
+tmp1 = permute(tensorprod(Eval,coeffs,2,1),[2 3 1 4]);
+tmp2 = permute(tensorprod(Eval,tmp1,2,1),[2 3 1 4]);
+vals = squeeze(permute(tensorprod(Eval,tmp2,2,1),[2 3 1 4]));
+
+end
\ No newline at end of file
diff --git a/@treefun3/coeffs2vals.m b/@treefun3/coeffs2vals.m
new file mode 100644
index 0000000..b399bee
--- /dev/null
+++ b/@treefun3/coeffs2vals.m
@@ -0,0 +1,23 @@
+function vals = coeffs2vals(coeffs)
+%
+%
+
+persistent Eval pstored
+[p,~,~,nd] = size(coeffs);
+
+if ( isempty(Eval) || p ~= pstored )
+    pstored = p;
+    Eval = ones(p, p);
+    x = linspace(-1, 1, p).';
+    x = (1-cos(pi*(2*(1:p)'-1)/(2*p)))/2;
+    Eval(:,2) = x;
+    for k=3:p
+        Eval(:,k) = 2*x.*Eval(:,k-1)-Eval(:,k-2);
+    end
+end
+
+tmp1 = permute(tensorprod(Eval,coeffs,2,1),[2 3 1 4]);
+tmp2 = permute(tensorprod(Eval,tmp1,2,1),[2 3 1 4]);
+vals = squeeze(permute(tensorprod(Eval,tmp2,2,1),[2 3 1 4]));
+
+end
\ No newline at end of file
diff --git a/@treefun3/leaves.m b/@treefun3/leaves.m
new file mode 100644
index 0000000..ae0cc4c
--- /dev/null
+++ b/@treefun3/leaves.m
@@ -0,0 +1,7 @@
+function ids = leaves(f)
+%LEAVES   Get the leaf IDs in a TREEFUN2.
+%   LEAVES(F) returns the leaf IDs in the TREEFUN2 F.
+
+ids = f.id(f.height == 0);
+
+end
diff --git a/@treefun3/plot.m b/@treefun3/plot.m
new file mode 100644
index 0000000..407c04f
--- /dev/null
+++ b/@treefun3/plot.m
@@ -0,0 +1,290 @@
+function varargout = plot(f, varargin)
+%PLOT   Plot a TREEFUN3.
+%   PLOT(F) gives a 3D plot of the TREEFUN3 F and shows the tree on
+%   which F is defined.
+%
+%   See also SURF, MESH.
+
+func = varargin{1};
+
+holdState = ishold();
+nplotpts = 51;
+
+h = instantiateSlice3GUI();
+handles = guihandles(h);
+
+% for k = 1:length(ids)
+%     id = ids(k);
+%     [x, y] = meshgrid(linspace(f.domain(1,id), f.domain(2,id), nplotpts), ...
+%                       linspace(f.domain(3,id), f.domain(4,id), nplotpts));
+%     u = coeffs2plotvals(f.coeffs{id});
+%     hk = surface(x, y, 0*u, u, 'EdgeAlpha', 1, varargin{:});
+%     if ( nargout > 0 )
+%         h(k) = hk; %#ok<AGROW>
+%     end
+% end
+
+[xx, yy, zz] = meshgrid(linspace(f.domain(1), f.domain(2), nplotpts), ...
+                        linspace(f.domain(3), f.domain(4), nplotpts), ...
+                        linspace(f.domain(5), f.domain(6), nplotpts));
+% v = func( xx, yy, zz);
+tmpval = func(0,0,0);
+nd = numel(tmpval);
+% nd = numel(func);
+% if nd == 1 && ~iscell(func)
+%   func1 = []; 
+%   func1{1} = func;
+%   func = func1;
+% end
+v = zeros(size(xx));
+tmpvals = func(xx,yy,zz);
+for k = 1:nd
+  v = v + tmpvals(:,:,:,k);
+end
+% for k = 1:nd
+%   v = v + func{k}(xx,yy,zz);
+% end
+if isreal(v)
+    [row,col,tube] = ind2sub(size(v), find(v(:) == max(v(:)), 1, 'last'));
+else
+    [row, col, tube] = ind2sub(size(v), find(abs(v(:)) == max(abs(v(:))), 1, 'last'));    
+end
+xslice = xx(row,col,tube); 
+yslice = yy(row,col,tube); 
+zslice = zz(row,col,tube); 
+
+set(handles.xSlider, 'Min', f.domain(1));
+set(handles.xSlider, 'Max', f.domain(2));
+set(handles.xSlider, 'Value', xslice);
+
+set(handles.ySlider, 'Min', f.domain(3));
+set(handles.ySlider, 'Max', f.domain(4));
+set(handles.ySlider, 'Value', yslice);
+
+set(handles.zSlider, 'Min', f.domain(5));
+set(handles.zSlider, 'Max', f.domain(6));
+set(handles.zSlider, 'Value', zslice);
+
+nSteps = 15; % number of slices allowed
+set(handles.xSlider, 'SliderStep', [1/nSteps , 1 ]);
+set(handles.ySlider, 'SliderStep', [1/nSteps , 1 ]);
+set(handles.zSlider, 'SliderStep', [1/nSteps , 1 ]);
+
+% Choose default command line output for the slice command:
+handles.xx = xx;
+handles.yy = yy;
+handles.zz = zz;
+handles.xslice = xslice;
+handles.yslice = yslice;
+handles.zslice = zslice;
+handles.v = v;
+
+if isreal(v)
+    slice(xx, yy, zz, v, xslice, yslice, zslice)
+    shading interp
+    colorbar
+else
+    hh = slice(xx, yy, zz, angle(-v), xslice, yslice, zslice); 
+    set(hh, 'EdgeColor','none')
+    caxis([-pi pi]),
+    colormap('hsv')
+    axis('equal')     
+end
+
+hold on
+
+% Plot the boxes
+ids = leaves(f);
+xdata = [f.domain([1 2 2 1 1 1 2 2 1 1], ids) ; nan(1, length(ids)); ... % bottom & top
+         f.domain([2 2 2 2 1 1], ids) ; nan(1, length(ids));]; % the rest to complete the box
+ydata = [f.domain([3 3 4 4 3 3 3 4 4 3], ids) ; nan(1, length(ids)); ...
+         f.domain([3 3 4 4 4 4], ids) ; nan(1, length(ids));];
+zdata = [f.domain([5 5 5 5 5 6 6 6 6 6], ids) ; nan(1, length(ids)); ...
+         f.domain([5 6 6 5 5 6], ids) ; nan(1, length(ids));];
+line('XData', xdata(:), 'YData', ydata(:), 'ZData', zdata(:), 'LineWidth', 1)
+hold off
+
+axis equal
+xlim(f.domain(1:2, 1))
+ylim(f.domain(3:4, 1))
+zlim(f.domain(5:6, 1))
+
+%
+handles.xdata = xdata(:);
+handles.ydata = ydata(:);
+handles.zdata = zdata(:);
+
+% Update handles structure
+guidata(h, handles);
+handles.output = handles.xSlider;
+
+if ( ~holdState )
+    hold off
+end
+
+if ( nargout > 0 )
+    varargout = {h};
+end
+
+% Force the figure to clear when another plot is drawn on it so that GUI
+% widgets don't linger.  (NB:  This property needs to be reset to 'add' every
+% time we change the plot using a slider; otherwise, the slider movement will
+% itself clear the figure, which is not what we want.)
+set(h, 'NextPlot', 'replacechildren');
+
+end
+
+function h = instantiateSlice3GUI()
+
+% Load up the GUI from the *.fig file.
+installDir = treefunroot();
+h = openFigInCurrentFigure([installDir '/@treefun3/slice.fig']);
+
+% Do any required initialization of the handle graphics objects.
+G = get(h, 'Children');
+for (i = 1:1:length(G))
+    if ( isa(G(i), 'matlab.ui.control.UIControl') )
+        % Adjust the background colors of the sliders.
+        if ( strcmp(G(i).Style, 'slider') )
+            if ( isequal(get(G(i), 'BackgroundColor'), ...
+                    get(0, 'defaultUicontrolBackgroundColor')) )
+                set(G(i), 'BackgroundColor', [.9 .9 .9]);
+            end
+        end
+        % Register callbacks.
+        switch ( G(i).Tag )
+            case 'xSlider'
+                G(i).Callback = @(hObj, data) ...
+                    xSlider_Callback(hObj, data, guidata(hObj));
+            case 'ySlider'
+                G(i).Callback = @(hObj, data) ...
+                    ySlider_Callback(hObj, data, guidata(hObj));
+            case 'zSlider'
+                G(i).Callback = @(hObj, data) ...
+                    zSlider_Callback(hObj, data, guidata(hObj));
+        end
+    end
+end
+
+% Store handles to GUI objects so that the callbacks can access them. 
+guidata(h, guihandles(h));
+
+end
+
+% --- Executes on xSlider movement.
+function xSlider_Callback(hObject, eventdata, handles)
+% hObject    handle to xSlider (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+nextPlot = get(hObject.Parent, 'NextPlot');
+set(hObject.Parent, 'NextPlot', 'add');
+
+xslice = get(hObject, 'Value');         %returns position of slider
+yslice = get(handles.ySlider, 'Value'); %returns position of slider
+zslice = get(handles.zSlider, 'Value'); %returns position of slider
+
+% The next slice command clears the title, if there was any. So, get that
+% and put it again afterwards.
+tit = get(gca(), 'Title');
+titText = tit.String;
+
+if ( isreal(handles.v) )
+    handles.slice = slice(handles.xx, handles.yy, handles.zz, handles.v, ...
+        xslice, yslice, zslice);
+    shading interp
+    colorbar, 
+else
+    handles.slice = slice(handles.xx, handles.yy, handles.zz, angle(-handles.v), ...
+        xslice, yslice, zslice); 
+    set(handles.slice, 'EdgeColor','none')
+    caxis([-pi pi]),
+    colormap('hsv')
+    axis('equal')
+end
+handles.line = line('XData', handles.xdata(:), 'YData', handles.ydata(:), 'ZData', handles.zdata(:), 'LineWidth', 1);
+title(titText)
+handles.output = hObject;
+
+set(hObject.Parent, 'NextPlot', nextPlot);
+
+end
+
+function ySlider_Callback(hObject, eventdata, handles)
+% --- Executes on ySlider movement.
+% hObject    handle to ySlider (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+nextPlot = get(hObject.Parent, 'NextPlot');
+set(hObject.Parent, 'NextPlot', 'add');
+
+yslice = get(hObject, 'Value');         %returns position of slider
+xslice = get(handles.xSlider, 'Value'); %returns position of slider
+zslice = get(handles.zSlider, 'Value'); %returns position of slider
+
+% The next slice command clears the title, if there was any. So, get that
+% and put it again afterwards.
+tit = get(gca(), 'Title');
+titText = tit.String;
+
+if ( isreal(handles.v) )
+    handles.slice = slice(handles.xx, handles.yy, handles.zz, handles.v, ...
+        xslice, yslice, zslice);
+    shading interp
+    colorbar, 
+else
+    handles.slice = slice(handles.xx, handles.yy, handles.zz, angle(-handles.v), ...
+        xslice, yslice, zslice); 
+    set(handles.slice, 'EdgeColor','none')
+    caxis([-pi pi]),
+    colormap('hsv')
+    axis('equal')    
+end
+handles.line = line('XData', handles.xdata(:), 'YData', handles.ydata(:), 'ZData', handles.zdata(:), 'LineWidth', 1);
+title(titText)
+handles.output = hObject;
+
+set(hObject.Parent, 'NextPlot', nextPlot);
+
+end
+
+function zSlider_Callback(hObject, eventdata, handles)
+% --- Executes on zSlider movement.
+% hObject    handle to zSlider (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+nextPlot = get(hObject.Parent, 'NextPlot');
+set(hObject.Parent, 'NextPlot', 'add');
+
+zslice = get(hObject, 'Value');         %returns position of slider
+xslice = get(handles.xSlider, 'Value'); %returns position of slider
+yslice = get(handles.ySlider, 'Value'); %returns position of slider
+
+% The next slice command clears the title, if there was any. So, get that
+% and put it again afterwards.
+tit = get(gca(), 'Title');
+titText = tit.String;
+
+if ( isreal(handles.v) )
+    handles.slice = slice(handles.xx, handles.yy, handles.zz, handles.v, ...
+        xslice, yslice, zslice);
+    shading interp
+    colorbar
+else
+    handles.slice = slice(handles.xx, handles.yy, handles.zz, angle(-handles.v), ...
+        xslice, yslice, zslice); 
+    set(handles.slice, 'EdgeColor','none')
+    caxis([-pi pi]),
+    colormap('hsv')
+    axis('equal')    
+end
+handles.line = line('XData', handles.xdata, 'YData', handles.ydata, 'ZData', handles.zdata, 'LineWidth', 1);
+title(titText)
+handles.output = hObject;
+
+set(hObject.Parent, 'NextPlot', nextPlot);
+
+end
+
diff --git a/@treefun3/private/coeffs2children.m b/@treefun3/private/coeffs2children.m
new file mode 100644
index 0000000..c0db47d
--- /dev/null
+++ b/@treefun3/private/coeffs2children.m
@@ -0,0 +1,89 @@
+function [LLD, LRD, ULD, URD, LLT, LRT, ULT, URT] = coeffs2children3d(coeffs)
+%COEFFS2CHILDREN   Convert 3D Chebyshev coefficients on a parent to 2D
+% Chebyshev coefficients on its four children.
+% follow Z curve, double check with ngrid
+
+persistent EvalL EvalR Fp pstored
+p = size(coeffs, 1);
+
+if ( isempty(EvalL) || isempty(EvalR) || isempty(Fp) || p ~= pstored )
+    pstored = p;
+    x = -cos(pi*(2*(1:p)'-1)/(2*p));
+    xL = 1/2*x-1/2;
+    xR = 1/2*x+1/2;
+    EvalL = ones(p, p); EvalR = ones(p, p); 
+    EvalL(:,2) = xL; EvalR(:,2) = xR;
+    for k=3:p
+      EvalL(:,k) = 2*xL.*EvalL(:,k-1)-EvalL(:,k-2);
+      EvalR(:,k) = 2*xR.*EvalR(:,k-1)-EvalR(:,k-2);
+    end
+    Fp = 2*cos(pi*((1:p)-1)'*(2*(p:-1:1)-1)/(2*p))/p;
+    Fp(1,:) = 1/2*Fp(1,:);
+end
+
+% Lower left down
+tmp1 = permute(tensorprod(EvalL,coeffs,2,1),[2 3 1 4]);
+tmp2 = permute(tensorprod(EvalL,tmp1,2,1),[2 3 1 4]);
+valsLLD = squeeze(permute(tensorprod(EvalL,tmp2,2,1),[2 3 1 4]));
+tmp1hat = permute(tensorprod(Fp,valsLLD,2,1),[2 3 1 4]);
+tmp2hat = permute(tensorprod(Fp,tmp1hat,2,1),[2 3 1 4]);
+LLD     = permute(tensorprod(Fp,tmp2hat,2,1),[2 3 1 4]);
+% Lower right down
+% [xx0LRD, yy0LRD, zz0LRD] = ndgrid(xR,xL,xL);
+tmp1 = permute(tensorprod(EvalR,coeffs,2,1),[2 3 1 4]);
+tmp2 = permute(tensorprod(EvalL,tmp1,2,1),[2 3 1 4]);
+valsLRD = squeeze(permute(tensorprod(EvalL,tmp2,2,1),[2 3 1 4]));
+tmp1hat = permute(tensorprod(Fp,valsLRD,2,1),[2 3 1 4]);
+tmp2hat = permute(tensorprod(Fp,tmp1hat,2,1),[2 3 1 4]);
+LRD     = permute(tensorprod(Fp,tmp2hat,2,1),[2 3 1 4]);
+% Upper left down
+% [xx0ULD, yy0ULD, zz0ULD] = ndgrid(xL,xR,xL);
+tmp1 = permute(tensorprod(EvalL,coeffs,2,1),[2 3 1 4]);
+tmp2 = permute(tensorprod(EvalR,tmp1,2,1),[2 3 1 4]);
+valsULD = squeeze(permute(tensorprod(EvalL,tmp2,2,1),[2 3 1 4]));
+tmp1hat = permute(tensorprod(Fp,valsULD,2,1),[2 3 1 4]);
+tmp2hat = permute(tensorprod(Fp,tmp1hat,2,1),[2 3 1 4]);
+ULD     = permute(tensorprod(Fp,tmp2hat,2,1),[2 3 1 4]);
+% Upper right down
+% [xx0URD, yy0URD, zz0URD] = ndgrid(xR,xR,xL);
+tmp1 = permute(tensorprod(EvalR,coeffs,2,1),[2 3 1 4]);
+tmp2 = permute(tensorprod(EvalR,tmp1,2,1),[2 3 1 4]);
+valsURD = squeeze(permute(tensorprod(EvalL,tmp2,2,1),[2 3 1 4]));
+tmp1hat = permute(tensorprod(Fp,valsURD,2,1),[2 3 1 4]);
+tmp2hat = permute(tensorprod(Fp,tmp1hat,2,1),[2 3 1 4]);
+URD     = permute(tensorprod(Fp,tmp2hat,2,1),[2 3 1 4]);
+
+% Lower left top
+% [xx0LLT, yy0LLT, zz0LLT] = ndgrid(xL,xL,xR);
+tmp1 = permute(tensorprod(EvalL,coeffs,2,1),[2 3 1 4]);
+tmp2 = permute(tensorprod(EvalL,tmp1,2,1),[2 3 1 4]);
+valsLLT = squeeze(permute(tensorprod(EvalR,tmp2,2,1),[2 3 1 4]));
+tmp1hat = permute(tensorprod(Fp,valsLLT,2,1),[2 3 1 4]);
+tmp2hat = permute(tensorprod(Fp,tmp1hat,2,1),[2 3 1 4]);
+LLT     = permute(tensorprod(Fp,tmp2hat,2,1),[2 3 1 4]);
+% Lower right top
+% [xx0LRT, yy0LRT, zz0LRT] = ndgrid(xR,xL,xR);
+tmp1 = permute(tensorprod(EvalR,coeffs,2,1),[2 3 1 4]);
+tmp2 = permute(tensorprod(EvalL,tmp1,2,1),[2 3 1 4]);
+valsLRT = squeeze(permute(tensorprod(EvalR,tmp2,2,1),[2 3 1 4]));
+tmp1hat = permute(tensorprod(Fp,valsLRT,2,1),[2 3 1 4]);
+tmp2hat = permute(tensorprod(Fp,tmp1hat,2,1),[2 3 1 4]);
+LRT     = permute(tensorprod(Fp,tmp2hat,2,1),[2 3 1 4]);
+% Upper left top
+% [xx0ULT, yy0ULT, zz0ULT] = ndgrid(xL,xR,xR);
+tmp1 = permute(tensorprod(EvalL,coeffs,2,1),[2 3 1 4]);
+tmp2 = permute(tensorprod(EvalR,tmp1,2,1),[2 3 1 4]);
+valsULT = squeeze(permute(tensorprod(EvalR,tmp2,2,1),[2 3 1 4]));
+tmp1hat = permute(tensorprod(Fp,valsULT,2,1),[2 3 1 4]);
+tmp2hat = permute(tensorprod(Fp,tmp1hat,2,1),[2 3 1 4]);
+ULT     = permute(tensorprod(Fp,tmp2hat,2,1),[2 3 1 4]);
+% Upper right top
+% [xx0URT, yy0URT, zz0URT] = ndgrid(xR,xR,xR);
+tmp1 = permute(tensorprod(EvalR,coeffs,2,1),[2 3 1 4]);
+tmp2 = permute(tensorprod(EvalR,tmp1,2,1),[2 3 1 4]);
+valsURT = squeeze(permute(tensorprod(EvalR,tmp2,2,1),[2 3 1 4]));
+tmp1hat = permute(tensorprod(Fp,valsURT,2,1),[2 3 1 4]);
+tmp2hat = permute(tensorprod(Fp,tmp1hat,2,1),[2 3 1 4]);
+URT     = permute(tensorprod(Fp,tmp2hat,2,1),[2 3 1 4]);
+
+end
diff --git a/@treefun3/private/isLeaf.m b/@treefun3/private/isLeaf.m
new file mode 100644
index 0000000..044e4c9
--- /dev/null
+++ b/@treefun3/private/isLeaf.m
@@ -0,0 +1,6 @@
+function out = isLeaf(f, id)
+%ISLEAF   Is this box a leaf?
+
+out = ( f.height(id) == 0 );
+
+end
diff --git a/@treefun3/private/openFigInCurrentFigure.m b/@treefun3/private/openFigInCurrentFigure.m
new file mode 100644
index 0000000..ff37808
--- /dev/null
+++ b/@treefun3/private/openFigInCurrentFigure.m
@@ -0,0 +1,35 @@
+function h = openFigInCurrentFigure(figfile)
+%OPENFIGINCURRENTFIGURE    Wrapper for OPENFIG to open in current figure.
+%   H = OPENFIGINCURRENTFIGURE(FIGFILE) calls OPENFIG to open the *.fig
+%   figure specified by the string FIGFILE but does some additional work
+%   to force it to be opened in the current figure instead of in a new one.
+
+% Copyright 2017 by The University of Oxford and The Chebfun Developers.
+% See http://www.chebfun.org/ for Chebfun information.
+
+% Developer note:  This function exists to allow a user to force the
+% special GUI-based plots for chebfun3 (e.g., isosurface(), slice(), etc.)
+% to appear in a figure of choice (e.g., by calling figure() to set the
+% current figure) instead of appearing in a new figure every time, just
+% like ordinary MATLAB plots.
+
+% Get rid of everything in the current figure.
+hgcf = gcf();
+clf();
+
+% Reset the 'NextPlot' property.  (Should we call newplot() here instead?)
+set(hgcf, 'NextPlot', 'add');
+
+% Open the *.fig file in a hidden window.
+h = openfig(figfile, 'invisible');
+
+% Copy the graphics elements from the *.fig figure into the current figure.
+copyobj(allchild(h), hgcf);
+
+% We don't need the *.fig figure anymore, so delete it.
+delete(h);
+
+% Return a handle to the current figure.
+h = hgcf;
+
+end
diff --git a/@treefun3/refineBox.m b/@treefun3/refineBox.m
new file mode 100644
index 0000000..7ac0e61
--- /dev/null
+++ b/@treefun3/refineBox.m
@@ -0,0 +1,247 @@
+function f = refineBox(f, id, func)
+
+persistent xx0 yy0 zz0 ww0 nstored
+nalias = f.n;
+nd = numel(func);
+if ( isempty(xx0) || f.n ~= nstored )
+    nstored = f.n;
+    x0 = (1-cos(pi*(2*(1:nalias)'-1)/(2*nalias)))/2;
+    [xx0, yy0, zz0] = ndgrid(x0);
+    l = floor(nalias/2)+1;
+    v = [2*exp(1i*pi*(0:nalias-l)/nalias)./(1-4*(0:nalias-l).^2)  zeros(1,l)];
+    w0 = real(ifft(v(1:nalias) + conj(v(nalias+1:-1:2))))'/2;
+    [wx0, wy0, wz0] = ndgrid(w0); 
+    ww0 = wx0.*wy0.*wz0;
+end
+
+% Split into eight child boxes
+dom = f.domain(:,id);
+xmid = mean(dom(1:2));
+ymid = mean(dom(3:4));
+zmid = mean(dom(5:6));
+
+cid1 = length(f.id)+1;
+f.domain(:,cid1)   = [dom(1) xmid dom(3) ymid dom(5) zmid];
+f.id(cid1)         = cid1;
+f.parent(cid1)     = id;
+f.children(:,cid1) = 0;
+f.level(cid1)      = f.level(id)+1;
+f.height(cid1)     = 0;
+cdom1              = f.domain(:,cid1);
+csclx1             = diff(cdom1(1:2));
+cscly1             = diff(cdom1(3:4));
+csclz1             = diff(cdom1(5:6));
+cxx1               = csclx1*xx0 + cdom1(1); 
+cyy1               = cscly1*yy0 + cdom1(3); 
+czz1               = csclz1*zz0 + cdom1(5); 
+cvals1             = cell(1,nd);
+for k = 1:nd
+  cvals1{k}        = func{k}(cxx1,cyy1,czz1);
+end
+cvals1             = cat(4,cvals1{:}); % additional for nd
+ccoeffs1           = treefun3.vals2coeffs(cvals1);
+f.coeffs{cid1}     = ccoeffs1(1:f.n,1:f.n,1:f.n,:); % to replace f.coeffs{cid1} = [];
+f.rint(:,cid1)     = squeeze(sqrt((csclx1*cscly1*csclz1)*sum(cvals1.^2.*ww0, [1 2 3])));
+f.vmax(:,cid1)     = squeeze(max(abs(cvals1),[],[1 2 3]));
+f.col(cid1)        = 2*f.col(id);
+f.row(cid1)        = 2*f.row(id);
+f.dep(cid1)        = 2*f.dep(id);
+% f.morton(cid1)     = cartesian2morton(f.col(cid1), f.row(cid1));
+
+cid2 = length(f.id)+1;
+f.domain(:,cid2)   = [xmid dom(2) dom(3) ymid dom(5) zmid];
+f.id(cid2)         = cid2;
+f.parent(cid2)     = id;
+f.children(:,cid2) = 0;
+f.level(cid2)      = f.level(id)+1;
+f.height(cid2)     = 0;
+cdom2              = f.domain(:,cid2);
+csclx2             = diff(cdom2(1:2));
+cscly2             = diff(cdom2(3:4));
+csclz2             = diff(cdom2(5:6));
+cxx2               = csclx2*xx0 + cdom2(1); 
+cyy2               = cscly2*yy0 + cdom2(3); 
+czz2               = csclz2*zz0 + cdom2(5); 
+cvals2             = cell(1,nd);
+for k = 1:nd
+  cvals2{k}        = func{k}(cxx2,cyy2,czz2);
+end
+cvals2             = cat(4,cvals2{:}); % additional for nd
+ccoeffs2           = treefun3.vals2coeffs(cvals2);
+f.coeffs{cid2}     = ccoeffs2(1:f.n,1:f.n,1:f.n,:); 
+f.rint(:,cid2)     = squeeze(sqrt((csclx2*cscly2*csclz2)*sum(cvals2.^2.*ww0, [1 2 3])));
+f.vmax(:,cid2)     = squeeze(max(abs(cvals2),[],[1 2 3]));
+f.col(cid2)        = 2*f.col(id) + 1;
+f.row(cid2)        = 2*f.row(id);
+f.dep(cid2)        = 2*f.dep(id);
+% f.morton(cid2)     = cartesian2morton(f.col(cid2), f.row(cid2));
+
+cid3 = length(f.id)+1;
+f.domain(:,cid3)   = [dom(1) xmid ymid dom(4) dom(5) zmid];
+f.id(cid3)         = cid3;
+f.parent(cid3)     = id;
+f.children(:,cid3) = 0;
+f.level(cid3)      = f.level(id)+1;
+f.height(cid3)     = 0;
+cdom3              = f.domain(:,cid3);
+csclx3             = diff(cdom3(1:2));
+cscly3             = diff(cdom3(3:4));
+csclz3             = diff(cdom3(5:6));
+cxx3               = csclx3*xx0 + cdom3(1); 
+cyy3               = cscly3*yy0 + cdom3(3); 
+czz3               = csclz3*zz0 + cdom3(5); 
+cvals3             = cell(1,nd);
+for k = 1:nd
+  cvals3{k}        = func{k}(cxx3,cyy3,czz3);
+end
+cvals3             = cat(4,cvals3{:}); % additional for nd
+ccoeffs3           = treefun3.vals2coeffs(cvals3);
+f.coeffs{cid3}     = ccoeffs3(1:f.n,1:f.n,1:f.n,:); 
+f.rint(:,cid3)     = squeeze(sqrt((csclx3*cscly3*csclz3)*sum(cvals3.^2.*ww0, [1 2 3])));
+f.vmax(:,cid3)     = squeeze(max(abs(cvals3),[],[1 2 3]));
+f.col(cid3)        = 2*f.col(id);
+f.row(cid3)        = 2*f.row(id) + 1;
+f.dep(cid3)        = 2*f.dep(id);
+% f.morton(cid3)     = cartesian2morton(f.col(cid3), f.row(cid3));
+
+cid4 = length(f.id)+1;
+f.domain(:,cid4)   = [xmid dom(2) ymid dom(4) dom(5) zmid];
+f.id(cid4)         = cid4;
+f.parent(cid4)     = id;
+f.children(:,cid4) = 0;
+f.level(cid4)      = f.level(id)+1;
+f.height(cid4)     = 0;
+cdom4              = f.domain(:,cid4);
+csclx4             = diff(cdom4(1:2));
+cscly4             = diff(cdom4(3:4));
+csclz4             = diff(cdom4(5:6));
+cxx4               = csclx4*xx0 + cdom4(1); 
+cyy4               = cscly4*yy0 + cdom4(3); 
+czz4               = csclz4*zz0 + cdom4(5); 
+cvals4             = cell(1,nd);
+for k = 1:nd
+  cvals4{k}        = func{k}(cxx4,cyy4,czz4);
+end
+cvals4             = cat(4,cvals4{:}); % additional for nd
+ccoeffs4           = treefun3.vals2coeffs(cvals4);
+f.coeffs{cid4}     = ccoeffs4(1:f.n,1:f.n,1:f.n,:); 
+f.rint(:,cid4)     = squeeze(sqrt((csclx4*cscly4*csclz4)*sum(cvals4.^2.*ww0, [1 2 3])));
+f.vmax(:,cid4)     = squeeze(max(abs(cvals4),[],[1 2 3]));
+f.col(cid4)        = 2*f.col(id) + 1;
+f.row(cid4)        = 2*f.row(id) + 1;
+f.dep(cid4)        = 2*f.dep(id);
+% f.morton(cid4)     = cartesian2morton(f.col(cid4), f.row(cid4));
+
+cid5 = length(f.id)+1;
+f.domain(:,cid5)   = [dom(1) xmid dom(3) ymid zmid dom(6)];
+f.id(cid5)         = cid5;
+f.parent(cid5)     = id;
+f.children(:,cid5) = 0;
+f.level(cid5)      = f.level(id)+1;
+f.height(cid5)     = 0;
+cdom5              = f.domain(:,cid5);
+csclx5             = diff(cdom5(1:2));
+cscly5             = diff(cdom5(3:4));
+csclz5             = diff(cdom5(5:6));
+cxx5               = csclx5*xx0 + cdom5(1); 
+cyy5               = cscly5*yy0 + cdom5(3); 
+czz5               = csclz5*zz0 + cdom5(5); 
+cvals5             = cell(1,nd);
+for k = 1:nd
+  cvals5{k}        = func{k}(cxx5,cyy5,czz5);
+end
+cvals5             = cat(4,cvals5{:}); % additional for nd
+ccoeffs5           = treefun3.vals2coeffs(cvals5);
+f.coeffs{cid5}     = ccoeffs5(1:f.n,1:f.n,1:f.n,:); 
+f.rint(:,cid5)     = squeeze(sqrt((csclx5*cscly5*csclz5)*sum(cvals5.^2.*ww0, [1 2 3])));
+f.vmax(:,cid5)     = squeeze(max(abs(cvals5),[],[1 2 3]));
+f.col(cid5)        = 2*f.col(id);
+f.row(cid5)        = 2*f.row(id);
+f.dep(cid5)        = 2*f.dep(id) + 1;
+
+cid6 = length(f.id)+1;
+f.domain(:,cid6)   = [xmid dom(2) dom(3) ymid zmid dom(6)];
+f.id(cid6)         = cid6;
+f.parent(cid6)     = id;
+f.children(:,cid6) = 0;
+f.level(cid6)      = f.level(id)+1;
+f.height(cid6)     = 0;
+cdom6              = f.domain(:,cid6);
+csclx6             = diff(cdom6(1:2));
+cscly6             = diff(cdom6(3:4));
+csclz6             = diff(cdom6(5:6));
+cxx6               = csclx6*xx0 + cdom6(1); 
+cyy6               = cscly6*yy0 + cdom6(3); 
+czz6               = csclz6*zz0 + cdom6(5); 
+cvals6             = cell(1,nd);
+for k = 1:nd
+  cvals6{k}        = func{k}(cxx6,cyy6,czz6);
+end
+cvals6             = cat(4,cvals6{:}); % additional for nd
+ccoeffs6           = treefun3.vals2coeffs(cvals6);
+f.coeffs{cid6}     = ccoeffs6(1:f.n,1:f.n,1:f.n,:); 
+f.rint(:,cid6)     = squeeze(sqrt((csclx6*cscly6*csclz6)*sum(cvals6.^2.*ww0, [1 2 3])));
+f.vmax(:,cid6)     = squeeze(max(abs(cvals6),[],[1 2 3]));
+f.col(cid6)        = 2*f.col(id) + 1;
+f.row(cid6)        = 2*f.row(id);
+f.dep(cid6)        = 2*f.dep(id) + 1;
+
+cid7 = length(f.id)+1;
+f.domain(:,cid7)   = [dom(1) xmid ymid dom(4) zmid dom(6)];
+f.id(cid7)         = cid7;
+f.parent(cid7)     = id;
+f.children(:,cid7) = 0;
+f.level(cid7)      = f.level(id)+1;
+f.height(cid7)     = 0;
+cdom7              = f.domain(:,cid7);
+csclx7             = diff(cdom7(1:2));
+cscly7             = diff(cdom7(3:4));
+csclz7             = diff(cdom7(5:6));
+cxx7               = csclx7*xx0 + cdom7(1); 
+cyy7               = cscly7*yy0 + cdom7(3); 
+czz7               = csclz7*zz0 + cdom7(5); 
+cvals7             = cell(1,nd);
+for k = 1:nd
+  cvals7{k}        = func{k}(cxx7,cyy7,czz7);
+end
+cvals7             = cat(4,cvals7{:}); % additional for nd
+ccoeffs7           = treefun3.vals2coeffs(cvals7);
+f.coeffs{cid7}     = ccoeffs7(1:f.n,1:f.n,1:f.n,:); 
+f.rint(:,cid7)     = squeeze(sqrt((csclx7*cscly7*csclz7)*sum(cvals7.^2.*ww0, [1 2 3])));
+f.vmax(:,cid7)     = squeeze(max(abs(cvals7),[],[1 2 3]));
+f.col(cid7)        = 2*f.col(id);
+f.row(cid7)        = 2*f.row(id) + 1;
+f.dep(cid7)        = 2*f.dep(id) + 1;
+
+cid8 = length(f.id)+1;
+f.domain(:,cid8)   = [xmid dom(2) ymid dom(4) zmid dom(6)];
+f.id(cid8)         = cid8;
+f.parent(cid8)     = id;
+f.children(:,cid8) = 0;
+f.level(cid8)      = f.level(id)+1;
+f.height(cid8)     = 0;
+cdom8              = f.domain(:,cid8);
+csclx8             = diff(cdom8(1:2));
+cscly8             = diff(cdom8(3:4));
+csclz8             = diff(cdom8(5:6));
+cxx8               = csclx8*xx0 + cdom8(1); 
+cyy8               = cscly8*yy0 + cdom8(3); 
+czz8               = csclz8*zz0 + cdom8(5); 
+cvals8             = cell(1,nd);
+for k = 1:nd
+  cvals8{k}        = func{k}(cxx8,cyy8,czz8);
+end
+cvals8             = cat(4,cvals8{:}); % additional for nd
+ccoeffs8           = treefun3.vals2coeffs(cvals8);
+f.coeffs{cid8}     = ccoeffs8(1:f.n,1:f.n,1:f.n,:); 
+f.rint(:,cid8)     = squeeze(sqrt((csclx8*cscly8*csclz8)*sum(cvals8.^2.*ww0, [1 2 3])));
+f.vmax(:,cid8)     = squeeze(max(abs(cvals8),[],[1 2 3]));
+f.col(cid8)        = 2*f.col(id) + 1;
+f.row(cid8)        = 2*f.row(id) + 1;
+f.dep(cid8)        = 2*f.dep(id) + 1;
+
+f.children(:,id) = [cid1 cid2 cid3 cid4 cid5 cid6 cid7 cid8];
+f.height(id) = 1;
+f.coeffs{id} = [];
+
+end
\ No newline at end of file
diff --git a/@treefun3/slice.fig b/@treefun3/slice.fig
new file mode 100644
index 0000000..dc97e51
Binary files /dev/null and b/@treefun3/slice.fig differ
diff --git a/@treefun3/treefun3.m b/@treefun3/treefun3.m
new file mode 100644
index 0000000..fe80afa
--- /dev/null
+++ b/@treefun3/treefun3.m
@@ -0,0 +1,725 @@
+classdef treefun3  %#ok<*PROP,*PROPLC>
+%TREEFUN3   Piecewise polynomial on an adaptive quadtree.
+%   TREEFUN3(F) constructs a TREEFUN3 object representing the function F on
+%   the domain [-1, 1] x [-1, 1] x [-1, 1]. F may be a function handle.
+%   A TREEFUN3 is constructed by recursively subdividing the domain until
+%   each piece is well approximated by a trivariate polynomial of degree
+%   (N-1) x (N-1). The default is N = 16. 
+%
+%   TREEFUN3(F, N) uses piecewise polynomials of degree (N-1) x (N-1).
+%
+%   TREEFUN3(F, [A B C D E F]) specifies a domain [A, B] x [C, D] x [E, F]
+%   on which the TREEFUN3 is defined.
+%
+%   TREEFUN3(F, [A B C D E F], N) specifies both a degree and a domain.
+%
+%   TREEFUN3(F, TF) creates a TREEFUN3 approximation to F using the tree
+%   structure from a previously-defined TREEFUN3. No adaptive construction
+%   takes place; the function F is simply initialized on the grid inherited
+%   from TF.
+%
+%   not all preperties are used ...
+
+    properties
+
+        n = 16
+        domain
+        level
+        height
+        id
+        parent
+        children
+        coeffs
+        col
+        row
+        dep
+        morton
+        flatNeighbors
+        leafNeighbors
+        rint0
+        rint
+        vmax
+
+    end
+
+    methods
+
+        function f = treefun3(varargin)
+
+            if ( nargin < 1 )
+                return
+            end
+
+            dom = [-1 1 -1 1 -1 1];
+            opts = struct();
+            opts.balance = false; % fortran way exists, morton way intial attempt in BBBBBBBoston
+            opts.neighbors = false;
+            opts.tol = 1e-12;
+            opts.checkpts = [];
+            opts.ifcoeffs = true; 
+            opts.ifstorecoeffs = true;
+            
+            if ( nargin == 2 )
+                if ( isa(varargin{2}, 'treefun3') ) % TREEFUN3(F, TF)
+                    % We were given the tree structure
+                    f = varargin{2};
+                    func = varargin{1};
+                    if ( isnumeric(func) && isscalar(func) )
+                        func = @(x,y) func + 0*x;
+                    elseif ( isa(func, 'chebfun3') )
+                        func = @(x,y) feval(func, x, y);
+                    end
+                    % We just need to fill in the leaf coefficients
+                    L = leaves(f);
+                    [xx0, yy0] = chebpts2(f.n, f.n, [0 1 0 1]);
+                    for id = L(:).'
+                        dom = f.domain(:,id);
+                        sclx = diff(dom(1:2));
+                        scly = diff(dom(3:4));
+                        xx = sclx*xx0 + dom(1);
+                        yy = scly*yy0 + dom(3);
+                        vals = func(xx,yy);
+                        f.coeffs{id} = treefun3.vals2coeffs(vals);
+                    end
+                    return
+                elseif ( isscalar(varargin{2}) ) % TREEFUN3(F, N)
+                    f.n = varargin{2};
+                else
+                    dom = varargin{2};           % TREEFUN3(F, [A B C D E F])
+                end
+            elseif ( nargin == 3 )
+                dom = varargin{2};               % TREEFUN3(F, [A B C D E F], N)
+                f.n = varargin{3};
+            elseif ( nargin == 4 )
+                dom = varargin{2};               % TREEFUN3(F, [A B C D E F], N, OPTS)
+                f.n = varargin{3};
+                opts1 = varargin{4};
+                if ( isfield(opts1, 'balance') )
+                    opts.balance = opts1.balance;
+                end
+                if ( isfield(opts1, 'neighbors') )
+                    opts.neighbors = opts1.neighbors;
+                end
+                if ( isfield(opts1, 'tol') )
+                    opts.tol = opts1.tol;
+                end
+                if ( isfield(opts1, 'checkpts') )
+                    opts.checkpts = opts1.checkpts;
+                end
+                if ( isfield(opts1, 'ifcoeffs') )
+                    opts.ifcoeffs = opts1.ifcoeffs;
+                end
+            elseif ( nargin == 9 )
+                % TREEFUN3(DOMAIN, LEVEL, HEIGHT, ID, PARENT, CHILDREN,
+                %   COEFFS, COL, ROW)
+                [f.domain, f.level, f.height, f.id, f.parent, ...
+                    f.children, f.coeffs, f.col, f.row] = deal(varargin{:});
+                f.morton = cartesian2morton(f.col, f.row);
+                f.n = size(f.coeffs{end}, 1);
+                f = balance(f);
+                [f.flatNeighbors, f.leafNeighbors] = generateNeighbors(f);
+                return
+            end
+
+            func = varargin{1};
+            if ( isnumeric(func) && isscalar(func) )
+                func = @(x,y) func + 0*x;
+            elseif ( isa(func, 'chebfun3') )
+                func = @(x,y) feval(func, x, y);
+            end
+
+            % initialize vals, rint, coeffs... wrap this up?
+            nalias = f.n;
+            tmpval = func(0,0,0);
+            nd = numel(tmpval);
+            % if nd == 1 && ~iscell(func)
+            %   func1 = []; 
+            %   func1{1} = func;
+            %   func = func1;
+            % end
+            x0 = (1-cos(pi*(2*(1:nalias)'-1)/(2*nalias)))/2;
+            [xx0, yy0, zz0] = ndgrid(x0);
+            l = floor(nalias/2)+1;
+            v = [2*exp(1i*pi*(0:nalias-l)/nalias)./(1-4*(0:nalias-l).^2)  zeros(1,l)];
+            w0 = real(ifft(v(1:nalias) + conj(v(nalias+1:-1:2))))'/2;
+            [wx0, wy0, wz0] = ndgrid(w0); 
+            sclx = diff(dom(1:2));
+            scly = diff(dom(3:4));
+            sclz = diff(dom(5:6));
+            xx = sclx*xx0 + dom(1); 
+            yy = scly*yy0 + dom(3); 
+            zz = sclz*zz0 + dom(5); 
+            ww0 = wx0.*wy0.*wz0;
+            vals = cell(1,nd);
+            tmpvals = func(xx,yy,zz);
+            for k = 1:nd
+              vals{k} = tmpvals(:,:,:,k);
+            end
+            % for k = 1:nd
+            %   vals{k} = func{k}(xx,yy,zz);
+            % end
+            vals = cat(4,vals{:});
+            coeffs = treefun3.vals2coeffs(vals);
+            rint = max(squeeze(sqrt((sclx*scly*sclz)*sum(vals.^2.*ww0, [1 2 3]))),1e-16); % initialze l2
+            
+            % 
+            f.domain(:,1)   = dom(:);
+            f.level(1)      = 0;
+            f.height(1)     = 0;
+            f.id(1)         = 1;
+            f.parent(1)     = 0;
+            f.children(:,1) = zeros(8, 1);
+            f.coeffs{1}     = [];
+            f.col           = uint64(0);
+            f.row           = uint64(0);
+            f.dep           = uint64(0); % is this the correct generalization?
+            f.coeffs{1}     = coeffs(1:f.n,1:f.n,1:f.n,:); 
+            f.rint0(:,1)    = zeros(nd,1); 
+            f.rint(:,1)     = rint;
+            f.vmax(:,1)     = squeeze(max(abs(vals),[],[1 2 3]));
+            
+            
+
+            % f = buildDepthFirst(f, func);
+            f = buildBreadthFirst(f, func, opts.tol, opts.checkpts, opts.ifcoeffs, opts.ifstorecoeffs);
+            f.morton = cartesian2morton(f.col, f.row, f.dep);
+
+            % Now do level restriction
+            if ( opts.balance )
+                f = balance(f);
+                % f = balancef(f);
+            else
+                % Do a cumulative sum in reverse to correct the heights
+                for k = length(f.id):-1:1
+                    if ( ~isLeaf(f, k) )
+                        f.height(k) = 1 + max(f.height(f.children(:,k)));
+                    end
+                end
+            end
+
+            if ( opts.neighbors )
+                [f.flatNeighbors, f.leafNeighbors] = generateNeighbors(f);
+            end
+
+        end
+
+    end
+
+    methods ( Access = private )
+
+        f = refineBox(f, id, func);
+
+        function f = buildBreadthFirst(f, func, tol, checkpts, ifcoeffs, ifstorecoeffs)
+
+            % initialization a 2nd tree that's the same as the 1st one
+            rint = f.rint(:,1);
+            [~,~,~,nd] = size(f.coeffs{1});
+            resolved = isResolved(f.coeffs{1}, f.domain(:,1), f.n, tol, f.vmax(:,1), func, checkpts, ifcoeffs, rint);
+            f.rint0(:,1) = rint;
+            if resolved
+              return
+            else
+              refinecnt = 1;
+              idl_start = 1; 
+              idl = 1; % current level box id
+            end
+            if ~ifstorecoeffs  % save memory, once isResolved done
+              f.coeffs{1} = [];
+            end
+
+            while(refinecnt)
+            
+              % process previous level unresolved boxes, idl(1), idl(2),....
+              domc    = zeros(6,8*refinecnt); % allocate space for all children info of this level's boxes
+              coeffsc = cell(1,8*refinecnt);
+              rintc   = zeros(nd,8*refinecnt);
+              rint0c  = zeros(nd,8*refinecnt);
+              vmaxc   = zeros(nd,8*refinecnt);
+              parentc = zeros(1,8*refinecnt);
+              idc     = zeros(1,8*refinecnt);
+              levelc  = zeros(1,8*refinecnt);
+              colc    = uint64(zeros(1,8*refinecnt));
+              rowc    = uint64(zeros(1,8*refinecnt));
+              depc    = uint64(zeros(1,8*refinecnt));
+              for k = 1:refinecnt % 1 by 1
+                id = idl(k);
+                [domck,coeffsck,rintck,vmaxck,colck,rowck,depck] = refineBoxv2(f.domain(:,id), f.n, func, f.col(id), f.row(id), f.dep(id)); % just refine the box...
+                f.children(:,id) = idl_start+8*(k-1)+(1:8); % start to know 
+                f.height(id)     = 1;
+                % 1 to 8
+                cidx             = 8*(k-1)+(1:8);
+                domc(:,cidx)     = domck;
+                coeffsc(cidx)    = coeffsck;
+                rintc(:,cidx)    = rintck;
+                vmaxc(:,cidx)    = vmaxck;
+                parentc(cidx)    = id*ones(1,8); % these children's parent id
+                idc(cidx)        = idl_start+8*(k-1)+(1:8); % children self id
+                levelc(cidx)     = (f.level(id)+1)*ones(1,8);
+                colc(cidx)       = colck;
+                rowc(cidx)       = rowck;
+                depc(cidx)       = depck;
+              end
+              childrenc  = zeros(8,8*refinecnt); % don't know yet, a bunch of 0s to be concatenate to the end of f.children, will know when next level
+              heightc    = zeros(1,8*refinecnt); % tmp leaf box
+              f.domain   = cat(2,f.domain,domc); % concatenate children info
+              f.coeffs   = cat(2,f.coeffs,coeffsc);
+              f.rint     = cat(2,f.rint,rintc);
+              f.rint0    = cat(2,f.rint0,rint0c);
+              f.vmax     = cat(2,f.vmax,vmaxc);
+              f.parent   = cat(2,f.parent,parentc);
+              f.id       = cat(2,f.id,idc);
+              f.level    = cat(2,f.level,levelc);
+              f.children = cat(2,f.children,childrenc);
+              f.height   = cat(2,f.height,heightc);
+              f.col      = cat(2,f.col,colc);
+              f.row      = cat(2,f.row,rowc);
+              f.dep      = cat(2,f.dep,depc);
+            
+              % update rint, - num of idl parents contribution + 8*num of idl children contribution
+              rint = sqrt(rint.^2 - sum(f.rint(:,idl).^2,2) + sum(f.rint(:,(idl_start+1):(idl_start+8*refinecnt)).^2,2)); 
+            
+              % see if newly created boxes need to be refined next
+              unresolvedl = true(1,8*length(idl)); % potentially unresolved
+              idl0 = (idl_start+1):(idl_start+8*length(idl)); % current level id, from start to num of idl refined
+              for k = 1:8*length(idl)
+                id = idl_start + k;
+                resolved = isResolved(f.coeffs{id}, f.domain(:,id), f.n, tol, f.vmax(:,id), func, checkpts, ifcoeffs, rint);
+                f.rint0(:,id) = rint;
+                if resolved
+                  unresolvedl(k) = false;
+                end
+                if ~ifstorecoeffs  % save memory, once isResolved done
+                  f.coeffs{id} = [];
+                end
+              end
+            
+              % for next loop
+              idl_start = idl_start + 8*refinecnt; % next level starts from here
+              idl = idl0(unresolvedl); % select subset of all potential boxes, based on unresolved
+              refinecnt = length(idl);
+              
+            end
+          
+            % % Note: the length changes at each iteration here
+            % id = 1;
+            % rint = f.rint(:,1);
+            % while ( id <= length(f.id) )
+            %     resolved = isResolved(f.coeffs{id}, f.domain(:,id), f.n, tol, f.vmax(:,id), func, checkpts, ifcoeffs, rint);
+            %     if ( resolved )
+            %         f.height(id) = 0;
+            %     else
+            %         % Split into eight child boxes
+            %         f = refineBox(f, id, func);
+            %         f.height(id) = 1;
+            %         f.coeffs{id} = []; % delete coeffs
+            %         rint = sqrt(rint.^2 - f.rint(:,id).^2 + sum(f.rint(:,end-7:end).^2,2)); % whenever split, update rint to be more accurate ...
+            %     end
+            %     id = id + 1;
+            % end
+
+            % Do a cumulative sum in reverse to correct the heights
+            for k = length(f.id):-1:1
+                if ( ~isLeaf(f, k) )
+                    %f.height(k) = f.height(k) + max(f.height(f.children(:,k)));
+                    f.height(k) = 1 + max(f.height(f.children(:,k)));
+                end
+            end
+
+        end
+
+    end
+
+    methods ( Static )
+
+        % u = poisson(f, isource);
+        coeffs = vals2coeffs(vals);
+        vals = coeffs2vals(coeffs);
+        vals = coeffs2refvals(coeffs);
+        % refvals = chebvals2refvals(chebvals);
+        checkvals = coeffs2checkvals(coeffs,x,y,z);
+
+    end
+    
+
+end
+
+function resolved = isResolved(coeffs, dom, n, tol, vmax, f, checkpts, ifcoeffs, rint)
+% need f if checkpts
+
+persistent xxx0 yyy0 zzz0 www0 nstored
+
+% nalias = n;
+nrefpts = n; % Sample at equispaced points to test error
+
+if ( isempty(xxx0) || n ~= nstored )
+    nstored = n;
+    [xxx0, yyy0, zzz0] = ndgrid(linspace(0, 1, nrefpts));
+    [wxxx0, wyyy0, wzzz0] = ndgrid(1/nstored*ones(nstored,1));
+    www0 = wxxx0.*wyyy0.*wzzz0; % weight
+end
+
+sclx = diff(dom(1:2));
+scly = diff(dom(3:4));
+sclz = diff(dom(5:6));
+
+h = sclx;
+eta = 0;
+[~,~,~,nd] = size(coeffs); % or rint...
+
+if ~ifcoeffs
+  resolved = 1;
+  if (sclx>tol) || (scly>tol) || (sclz>tol) % is this ok?
+    xxx = sclx*xxx0 + dom(1); 
+    yyy = scly*yyy0 + dom(3); 
+    zzz = sclz*zzz0 + dom(5); 
+    F = f(xxx,yyy,zzz);
+    G = treefun3.coeffs2refvals(coeffs); % structured grid
+    GmF = (G - F).^2 .* www0;
+    erra = sqrt(sum(reshape(GmF,[],nd),1));
+    resolved = prod(erra(:) < tol * sqrt(1/(sclx*scly*sclz)) * rint(:));
+
+    % vals = cell(1,nd);
+    % tmpvals = f(xxx,yyy,zzz);
+    % for k = 1:nd
+    %   vals{k} = tmpvals(:,:,:,k);
+    % end
+    % % for k = 1:nd
+    % %   vals{k} = f{k}(xxx,yyy,zzz);
+    % % end
+    % F = cat(4,vals{:});
+    % G = treefun3.coeffs2refvals(coeffs); % structured grid
+    % for k = 1:nd
+    %   erra = sqrt(sum(squeeze(G(:,:,:,k) - F(:,:,:,k)).^2.*www0,'all'));
+    %   resolved = resolved && ( erra < tol * sqrt(1/(sclx*scly*sclz)) * rint(k) );
+    % end
+    
+  end
+
+elseif ifcoeffs
+  resolved = 1;
+  erra = zeros(nd,1);
+  for k = 1:nd
+    erra(k) = sqrt(  sum(coeffs(end,:,:,k).^2,'all') ...
+                + sum(coeffs(1:end-1,end,:,k).^2,'all') ...
+                + sum(coeffs(1:end-1,1:end-1,end,k).^2,'all')) / sqrt(n^3 - (n-1)^3); % only last slice
+    % erra = sqrt(  sum(coeffs(end-1:end,:,:,k).^2,'all') ...
+    %             + sum(coeffs(1:end-2,end-1:end,:,k).^2,'all') ...
+    %             + sum(coeffs(1:end-2,1:end-2,end-1:end,k).^2,'all')) / sqrt(n^3 - (n-2)^3); 
+    resolved = resolved && ( erra(k) < tol* sqrt(1/(sclx*scly*sclz)) * rint(k) );
+  end
+
+else % did not check
+  vmax = max(abs(vals(:))); % if needed, compute and store when refineBox
+
+  Ex = sum(abs(coeffs(end-1:end,:,:)), 'all') / (3*n^2);
+  Ey = sum(abs(coeffs(:,end-1:end,:)), 'all') / (3*n^2);
+  Ez = sum(abs(coeffs(:,:,end-1:end)), 'all') / (3*n^2);
+  err_cfs = (Ex + Ey + Ez) / 3;
+  
+  % F = f(xxx,yyy,zzz);
+  % G = coeffs2refvals(coeffs);
+  % err_vals = max(abs(F(:) - G(:)));
+  
+  err = err_cfs;
+  %err = min(err_cfs, err_vals);
+  
+  resolved = ( err * h^eta < tol * max(vmax, 1) );
+end
+
+if ( ~isempty(checkpts) ) % check if func values @ checkpts agree
+  % func = @(x,y,z) exp(-(x.^2+y.^2+z.^2)*5000000) + exp(-((x-1/2).^2+(y-1/3).^2+(z-3/5).^2)*1000000) + exp(-((x+1/2).^2+(y+1/3).^2+(z+3/5).^2)*2000000);
+  % f = treefun3(func,[-2 2 -2 2 -2 2],10,struct('checkpts',[0 1/2 -1/2;0 1/3 -1/3;0 3/5 -3/5])); vs f = treefun3(func,[-2 2 -2 2 -2 2],10);
+  % plot(f,func)
+  % for later, reduce checkpts if resolved
+  xxx = 2 * ((checkpts(1,:)' - dom(1))/sclx) - 1;
+  yyy = 2 * ((checkpts(2,:)' - dom(3))/scly) - 1;
+  zzz = 2 * ((checkpts(3,:)' - dom(5))/sclz) - 1;
+  in = ( xxx>=-1 & xxx<=1 & ...
+         yyy>=-1 & yyy<=1 & ...
+         zzz>=-1 & zzz<=1);
+  if ( any(in) )
+    F = zeros(nd,sum(in));
+    F = reshape(squeeze(f(checkpts(1,in)',checkpts(2,in)',checkpts(3,in)'))',[nd sum(in)]);
+    % for k = 1:nd
+    %   F(k,:) = f{k}(checkpts(1,in)',checkpts(2,in)',checkpts(3,in)')'; % nd x n checkpts
+    % end
+    G = treefun3.coeffs2checkvals(coeffs,xxx(in),yyy(in),zzz(in));
+    err_checkvals = max(abs(F - G),[],2);
+    for k = 1:nd
+      resolved = resolved && ( err_checkvals(k) * h^eta < tol * max(vmax(k), 1) );
+    end
+  end
+end
+
+end
+
+
+function [domain,coeffs,rint,vmax,ccol,crow,cdep] = refineBoxv2(dom, n, func, col, row, dep)
+
+% if nargin==0, test_refineBox3dv2; return; end
+
+persistent xx0 yy0 zz0 ww0 nstored
+nalias = n;
+if ( isempty(xx0) || n ~= nstored )
+    nstored = n;
+    x0 = (1-cos(pi*(2*(1:nalias)'-1)/(2*nalias)))/2;
+    [xx0, yy0, zz0] = ndgrid(x0);
+    l = floor(nalias/2)+1;
+    v = [2*exp(1i*pi*(0:nalias-l)/nalias)./(1-4*(0:nalias-l).^2)  zeros(1,l)];
+    w0 = real(ifft(v(1:nalias) + conj(v(nalias+1:-1:2))))'/2;
+    [wx0, wy0, wz0] = ndgrid(w0); 
+    ww0 = wx0.*wy0.*wz0;
+end
+
+% Split into eight child boxes
+xmid = mean(dom(1:2));
+ymid = mean(dom(3:4));
+zmid = mean(dom(5:6));
+
+% domain and coeffs
+domain = zeros(6,8);
+coeffs = cell(1,8);
+% nd = numel(func);
+tmpval = func(0,0,0);
+nd = numel(tmpval);
+
+% morton related
+ccol = uint64(zeros(1,8));
+crow = uint64(zeros(1,8));
+cdep = uint64(zeros(1,8));
+
+
+cdom1              = [dom(1) xmid dom(3) ymid dom(5) zmid];
+csclx1             = diff(cdom1(1:2));
+cscly1             = diff(cdom1(3:4));
+csclz1             = diff(cdom1(5:6));
+cxx1               = csclx1*xx0 + cdom1(1); 
+cyy1               = cscly1*yy0 + cdom1(3); 
+czz1               = csclz1*zz0 + cdom1(5); 
+% cvals1             = cell(1,nd);
+% tmpvals1           = func(cxx1,cyy1,czz1);
+% for k = 1:nd
+%   cvals1{k}        = tmpvals1(:,:,:,k);
+% end
+% % for k = 1:nd
+% %   cvals1{k}        = func{k}(cxx1,cyy1,czz1);
+% % end
+% cvals1             = cat(4,cvals1{:});
+cvals1             = func(cxx1,cyy1,czz1);
+ccoeffs1           = treefun3.vals2coeffs(cvals1);
+ccol(1)            = 2*col;
+crow(1)            = 2*row;
+cdep(1)            = 2*dep;
+% f.morton(cid1)     = cartesian2morton(f.col(cid1), f.row(cid1));
+domain(:,1)        = cdom1;
+coeffs{1}          = ccoeffs1(1:n,1:n,1:n,:); % to replace f.coeffs{cid1} = [];
+
+cdom2              = [xmid dom(2) dom(3) ymid dom(5) zmid];
+csclx2             = diff(cdom2(1:2));
+cscly2             = diff(cdom2(3:4));
+csclz2             = diff(cdom2(5:6));
+cxx2               = csclx2*xx0 + cdom2(1); 
+cyy2               = cscly2*yy0 + cdom2(3); 
+czz2               = csclz2*zz0 + cdom2(5); 
+% cvals2             = cell(1,nd);
+% tmpvals2           = func(cxx2,cyy2,czz2);
+% for k = 1:nd
+%   cvals2{k}        = tmpvals2(:,:,:,k);
+% end
+% % for k = 1:nd
+% %   cvals2{k}        = func{k}(cxx2,cyy2,czz2);
+% % end
+% cvals2             = cat(4,cvals2{:});
+cvals2             = func(cxx2,cyy2,czz2);
+ccoeffs2           = treefun3.vals2coeffs(cvals2);
+ccol(2)            = 2*col + 1;
+crow(2)            = 2*row;
+cdep(2)            = 2*dep;
+% f.morton(cid2)     = cartesian2morton(f.col(cid2), f.row(cid2)); 
+domain(:,2)        = cdom2;
+coeffs{2}          = ccoeffs2(1:n,1:n,1:n,:); 
+
+cdom3              = [dom(1) xmid ymid dom(4) dom(5) zmid];
+csclx3             = diff(cdom3(1:2));
+cscly3             = diff(cdom3(3:4));
+csclz3             = diff(cdom3(5:6));
+cxx3               = csclx3*xx0 + cdom3(1); 
+cyy3               = cscly3*yy0 + cdom3(3); 
+czz3               = csclz3*zz0 + cdom3(5); 
+% cvals3             = cell(1,nd);
+% tmpvals3           = func(cxx3,cyy3,czz3);
+% for k = 1:nd
+%   cvals3{k}        = tmpvals3(:,:,:,k);
+% end
+% % for k = 1:nd
+% %   cvals3{k}        = func{k}(cxx3,cyy3,czz3);
+% % end
+% cvals3             = cat(4,cvals3{:});
+cvals3             = func(cxx3,cyy3,czz3);
+ccoeffs3           = treefun3.vals2coeffs(cvals3);
+ccol(3)            = 2*col;
+crow(3)            = 2*row + 1;
+cdep(3)            = 2*dep;
+% f.morton(cid3)     = cartesian2morton(f.col(cid3), f.row(cid3));
+domain(:,3)        = cdom3;
+coeffs{3}          = ccoeffs3(1:n,1:n,1:n,:); 
+
+cdom4              = [xmid dom(2) ymid dom(4) dom(5) zmid];
+csclx4             = diff(cdom4(1:2));
+cscly4             = diff(cdom4(3:4));
+csclz4             = diff(cdom4(5:6));
+cxx4               = csclx4*xx0 + cdom4(1); 
+cyy4               = cscly4*yy0 + cdom4(3); 
+czz4               = csclz4*zz0 + cdom4(5); 
+% cvals4             = cell(1,nd);
+% tmpvals4           = func(cxx4,cyy4,czz4);
+% for k = 1:nd
+%   cvals4{k}        = tmpvals4(:,:,:,k);
+% end
+% % for k = 1:nd
+% %   cvals4{k}        = func{k}(cxx4,cyy4,czz4);
+% % end
+% cvals4             = cat(4,cvals4{:});
+cvals4             = func(cxx4,cyy4,czz4);
+ccoeffs4           = treefun3.vals2coeffs(cvals4);
+ccol(4)            = 2*col + 1;
+crow(4)            = 2*row + 1;
+cdep(4)            = 2*dep;
+% f.morton(cid4)     = cartesian2morton(f.col(cid4), f.row(cid4));
+domain(:,4)        = cdom4;
+coeffs{4}          = ccoeffs4(1:n,1:n,1:n,:); 
+
+cdom5              = [dom(1) xmid dom(3) ymid zmid dom(6)];
+csclx5             = diff(cdom5(1:2));
+cscly5             = diff(cdom5(3:4));
+csclz5             = diff(cdom5(5:6));
+cxx5               = csclx5*xx0 + cdom5(1); 
+cyy5               = cscly5*yy0 + cdom5(3); 
+czz5               = csclz5*zz0 + cdom5(5); 
+% cvals5             = cell(1,nd);
+% tmpvals5           = func(cxx5,cyy5,czz5);
+% for k = 1:nd
+%   cvals5{k}        = tmpvals5(:,:,:,k);
+% end
+% % for k = 1:nd
+% %   cvals5{k}        = func{k}(cxx5,cyy5,czz5);
+% % end
+% cvals5             = cat(4,cvals5{:});
+cvals5             = func(cxx5,cyy5,czz5);
+ccoeffs5           = treefun3.vals2coeffs(cvals5);
+ccol(5)            = 2*col;
+crow(5)            = 2*row;
+cdep(5)            = 2*dep + 1;
+%        
+domain(:,5)        = cdom5;
+coeffs{5}          = ccoeffs5(1:n,1:n,1:n,:); 
+
+cdom6              = [xmid dom(2) dom(3) ymid zmid dom(6)];
+csclx6             = diff(cdom6(1:2));
+cscly6             = diff(cdom6(3:4));
+csclz6             = diff(cdom6(5:6));
+cxx6               = csclx6*xx0 + cdom6(1); 
+cyy6               = cscly6*yy0 + cdom6(3); 
+czz6               = csclz6*zz0 + cdom6(5); 
+% cvals6             = cell(1,nd);
+% tmpvals6           = func(cxx6,cyy6,czz6);
+% for k = 1:nd
+%   cvals6{k}        = tmpvals6(:,:,:,k);
+% end
+% % for k = 1:nd
+% %   cvals6{k}        = func{k}(cxx6,cyy6,czz6);
+% % end
+% cvals6             = cat(4,cvals6{:});
+cvals6             = func(cxx6,cyy6,czz6);
+ccoeffs6           = treefun3.vals2coeffs(cvals6);
+ccol(6)            = 2*col + 1;
+crow(6)            = 2*row;
+cdep(6)            = 2*dep + 1;
+%
+domain(:,6)        = cdom6;
+coeffs{6}          = ccoeffs6(1:n,1:n,1:n,:); 
+
+cdom7              = [dom(1) xmid ymid dom(4) zmid dom(6)];
+csclx7             = diff(cdom7(1:2));
+cscly7             = diff(cdom7(3:4));
+csclz7             = diff(cdom7(5:6));
+cxx7               = csclx7*xx0 + cdom7(1); 
+cyy7               = cscly7*yy0 + cdom7(3); 
+czz7               = csclz7*zz0 + cdom7(5); 
+% cvals7             = cell(1,nd);
+% tmpvals7           = func(cxx7,cyy7,czz7);
+% for k = 1:nd
+%   cvals7{k}        = tmpvals7(:,:,:,k);
+% end
+% % for k = 1:nd
+% %   cvals7{k}        = func{k}(cxx7,cyy7,czz7);
+% % end
+% cvals7             = cat(4,cvals7{:});
+cvals7             = func(cxx7,cyy7,czz7);
+ccoeffs7           = treefun3.vals2coeffs(cvals7);
+ccol(7)            = 2*col;
+crow(7)            = 2*row + 1;
+cdep(7)            = 2*dep + 1;
+%        
+domain(:,7)        = cdom7;
+coeffs{7}          = ccoeffs7(1:n,1:n,1:n,:); 
+
+cdom8              = [xmid dom(2) ymid dom(4) zmid dom(6)];
+csclx8             = diff(cdom8(1:2));
+cscly8             = diff(cdom8(3:4));
+csclz8             = diff(cdom8(5:6));
+cxx8               = csclx8*xx0 + cdom8(1); 
+cyy8               = cscly8*yy0 + cdom8(3); 
+czz8               = csclz8*zz0 + cdom8(5); 
+% cvals8             = cell(1,nd);
+% tmpvals8           = func(cxx8,cyy8,czz8);
+% for k = 1:nd
+%   cvals8{k}        = tmpvals8(:,:,:,k);
+% end
+% % for k = 1:nd
+% %   cvals8{k}        = func{k}(cxx8,cyy8,czz8);
+% % end
+% cvals8             = cat(4,cvals8{:});
+cvals8             = func(cxx8,cyy8,czz8);
+ccoeffs8           = treefun3.vals2coeffs(cvals8);
+ccol(8)            = 2*col + 1;
+crow(8)            = 2*row + 1;
+cdep(8)            = 2*dep + 1;
+domain(:,8)        = cdom8;
+coeffs{8}          = ccoeffs8(1:n,1:n,1:n,:); 
+
+% a few other things
+% rint = [squeeze(sqrt((csclx1*cscly1*csclz1)*sum(cvals1.^2.*ww0, [1 2 3]))),...
+%         squeeze(sqrt((csclx2*cscly2*csclz2)*sum(cvals2.^2.*ww0, [1 2 3]))),...
+%         squeeze(sqrt((csclx3*cscly3*csclz3)*sum(cvals3.^2.*ww0, [1 2 3]))),...
+%         squeeze(sqrt((csclx4*cscly4*csclz4)*sum(cvals4.^2.*ww0, [1 2 3]))),...
+%         squeeze(sqrt((csclx5*cscly5*csclz5)*sum(cvals5.^2.*ww0, [1 2 3]))),...
+%         squeeze(sqrt((csclx6*cscly6*csclz6)*sum(cvals6.^2.*ww0, [1 2 3]))),...
+%         squeeze(sqrt((csclx7*cscly7*csclz7)*sum(cvals7.^2.*ww0, [1 2 3]))),...
+%         squeeze(sqrt((csclx8*cscly8*csclz8)*sum(cvals8.^2.*ww0, [1 2 3])))];
+% vmax = [squeeze(max(abs(cvals1),[],[1 2 3])),...
+%         squeeze(max(abs(cvals2),[],[1 2 3])),...
+%         squeeze(max(abs(cvals3),[],[1 2 3])),...
+%         squeeze(max(abs(cvals4),[],[1 2 3])),...
+%         squeeze(max(abs(cvals5),[],[1 2 3])),...
+%         squeeze(max(abs(cvals6),[],[1 2 3])),...
+%         squeeze(max(abs(cvals7),[],[1 2 3])),...
+%         squeeze(max(abs(cvals8),[],[1 2 3]))];
+
+cvals = {cvals1, cvals2, cvals3, cvals4, cvals5, cvals6, cvals7, cvals8};
+csclx = [csclx1, csclx2, csclx3, csclx4, csclx5, csclx6, csclx7, csclx8];
+cscly = [cscly1, cscly2, cscly3, cscly4, cscly5, cscly6, cscly7, cscly8];
+csclz = [csclz1, csclz2, csclz3, csclz4, csclz5, csclz6, csclz7, csclz8];
+nsub = numel(cvals);
+rint = zeros(nd,nsub);
+vmax = zeros(nd,nsub);
+for i = 1:nsub
+  sw = csclx(i) * cscly(i) * csclz(i);
+  cvalstmp = reshape(cvals{i},[],nd);
+  rint(:,i) = (sw * sum(cvalstmp.*cvalstmp.*ww0(:),1)); 
+  vmax(:,i) = max(abs(cvalstmp)); 
+end
+rint = sqrt(rint);
+
+end
+
diff --git a/@treefun3/vals2coeffs.m b/@treefun3/vals2coeffs.m
new file mode 100644
index 0000000..a3b8e60
--- /dev/null
+++ b/@treefun3/vals2coeffs.m
@@ -0,0 +1,34 @@
+function coeffs = vals2coeffs(vals)
+%VALS2COEFFS   Convert values at Chebyshev points to Chebyshev
+% coefficients.
+%
+%   See also COEFFS2VALS.
+
+% Store the vals2coeffs matrices for sizes < cutoff
+persistent F
+cutoff = 1e+03; % ... later
+
+% Get the length of the input:
+[p,~,~,nd] = size(vals); % vals of dim p^3 x nd
+
+if ( p <= 1 )
+    % Trivial case (constant):
+    coeffs = vals;
+elseif ( p < cutoff )
+    % Use matrix multiplication for small problems
+    if ( isempty(F) )
+        F = cell(cutoff, 1);
+    end
+    if ( isempty(F{p}) )
+        F{p} = 2*cos(pi*((1:p)-1)'*(2*(p:-1:1)-1)/(2*p))/p;
+        F{p}(1,:) = 1/2*F{p}(1,:);
+    end
+    tmp1hat = permute(tensorprod(F{p},vals,2,1),[2 3 1 4]);
+    tmp2hat = permute(tensorprod(F{p},tmp1hat,2,1),[2 3 1 4]);
+    coeffs  = permute(tensorprod(F{p},tmp2hat,2,1),[2 3 1 4]);
+else
+    % Use fast transform ... not yet
+    % coeffs = chebtech2.vals2coeffs( chebtech2.vals2coeffs(vals).' ).';
+end
+
+end
diff --git a/cartesian2morton.m b/cartesian2morton.m
index fad4e4f..3959abe 100644
--- a/cartesian2morton.m
+++ b/cartesian2morton.m
@@ -1,6 +1,12 @@
-function morton = cartesian2morton(x, y)
+function morton = cartesian2morton(x, y, z)
 
-morton = bitshift(Part1By1_64(y), 1) + Part1By1_64(x);
+if nargin == 2
+  morton = bitshift(Part1By1_64(y), 1) + Part1By1_64(x);
+end
+
+if nargin == 3
+  morton = bitshift(Part1By2_64(z), 2) + bitshift(Part1By2_64(y), 1) + Part1By2_64(x);
+end
 
 end
 
@@ -15,3 +21,18 @@
 x = bitand(bitxor(x, bitshift(x, 1)),  0x5555555555555555);
 
 end
+
+function x = Part1By2_64(x)
+% "Insert" two 0 bits after each of the 16 low bits of x
+% https://github.com/trevorprater/pymorton/blob/f248683c19d90e904193c58bbd03e77bc2c43768/pymorton/pymorton.py#L20
+
+x = bitand(x, uint64(0x1fffff));
+x = bitand(bitxor(x, bitshift(x, 32)), uint64(0x1f00000000ffff));
+x = bitand(bitxor(x, bitshift(x, 16)), uint64(0x1f0000ff0000ff));
+x = bitand(bitxor(x, bitshift(x, 8)),  uint64(0x100f00f00f00f00f));
+x = bitand(bitxor(x, bitshift(x, 4)),  uint64(0x10c30c30c30c30c3));
+x = bitand(bitxor(x, bitshift(x, 2)),  uint64(0x1249249249249249));
+
+end
+
+