diff --git a/classes/heat.m b/classes/heat.m
index 3e1a2bb..36f39d9 100644
--- a/classes/heat.m
+++ b/classes/heat.m
@@ -44,7 +44,7 @@
             obj.heatDiffusion   = p.Results.heatDiffusion;
             obj.intpAtInterface = p.Results.intpAtInterface;
             % set default ode options after initialization of parent class
-            obj.odeOptions.RelTol = 1e-3;
+            obj.odeOptions.RelTol = 1e-4;
         end%function
         
         %% Display
@@ -635,7 +635,8 @@ function setBoundaryCondition(obj,varargin)
             index = finderb(z,dStart);
             unitCell = handles{index}; % this is the handle to the corresponding unitCell
             
-            k = cellfun(@feval,(unitCell.thermCond)',num2cell(T));            
+            k = cellfun(@feval,(unitCell.thermCond)',repmat({T},K,1));  % NTM
+%             k = cellfun(@feval,(unitCell.thermCond)',num2cell(T));      % 1TM
             % these are the parameters of the differential equation as they
             % are defined in matlab for the pdesolver
             
@@ -699,7 +700,7 @@ function setBoundaryCondition(obj,varargin)
         %% getEnergyMap
         % Returns an energy profile for the sample structure after
         % optical excitation.
-        function energyMap = getEnergyMap(obj,time,excitation,initTemp,tempMap)
+        function energyMap = getEnergyMap(obj,time,excitation,initTemp,tempMap,includeSubstrate)
             % create a unique hash
             hash = obj.getHash(time,excitation,initTemp);
             % create the file name to look for
@@ -710,7 +711,7 @@ function setBoundaryCondition(obj,varargin)
                 obj.dispMessage(['_energyMap_ loaded from file ' filename]);
             else
                 % file does not exist so calculate and save
-                energyMap = obj.calcEnergyMap(tempMap,initTemp);
+                energyMap = obj.calcEnergyMap(tempMap,initTemp,includeSubstrate);
                 save(filename,'energyMap');
                 obj.dispMessage(['_energyMap_ saved to file ' filename]);
             end                
@@ -719,7 +720,7 @@ function setBoundaryCondition(obj,varargin)
         %% calcEnergyMap
         % Calculates the energy profile for the sample structure after
         % optical excitation.
-        function energyMap = calcEnergyMap(obj,tempMap,initTemp)
+        function energyMap = calcEnergyMap(obj,tempMap,initTemp,includeSubstrate)
             
             disp('Calculating _energyMap_ ...')
             tic
@@ -731,10 +732,14 @@ function setBoundaryCondition(obj,varargin)
             aAxes           = obj.S.getUnitCellPropertyVector('aAxis');
             bAxes           = obj.S.getUnitCellPropertyVector('bAxis');
             UCmasses        = normMasses.*(aAxes/1e-10).*(bAxes/1e-10);             % calculates vector of unit cell masses
-            Cells           = obj.S.getNumberOfUnitCells;
+            if includeSubstrate
+                Cells           = obj.S.getNumberOfUnitCells;
+            else
+                Cells           = obj.S.getNumberOfUnitCells - length(obj.S.getAllPositionsPerUniqueUnitCell{end});
+            end
 
             for k=1:obj.S.numSubSystems
-                parfor i=1:size(tempMap,1)
+                for i=1:size(tempMap,1)
                     for n=1:Cells
                         energyMap(i,n,k) = UCmasses(n)*( intHeatCapacity{n,k}(tempMap(i,n,k)) - intHeatCapacity{n,k}(initTemp(n,k)) );
                     end
diff --git a/classes/phonon.m b/classes/phonon.m
index 62e0f37..fc36206 100644
--- a/classes/phonon.m
+++ b/classes/phonon.m
@@ -183,7 +183,10 @@ function disp(obj)
                     % current temperature for each subsystem                        
                     % traverse subsystems
                     for j=1:K
-                        intAlphaT(:,j) = cellfun(@feval,intLinThermExps(:,j),num2cell(squeeze(tempMap(i,:,j))'));
+                        % traverse unit cells
+                        for n=1:N
+                            intAlphaT(n,j) = intLinThermExps{n,j}(tempMap(i,n,j));      % cellfun(@feval,intLinThermExps(:,j),num2cell(squeeze(tempMap(i,:,j))'));
+                        end
                     end%for
 
                     % calculate the length of the sticks of each subsystem and sum
diff --git a/classes/unitCell.m b/classes/unitCell.m
index 7f53c85..4b0e0b8 100644
--- a/classes/unitCell.m
+++ b/classes/unitCell.m
@@ -111,6 +111,8 @@
                     obj.debWalFac{k} = str2func('@(T)(0.*T./T)');
                     obj.debWalFac(2) = debWalFacInit;
                 end%for
+            else
+                obj.debWalFac = debWalFacInit;
             end%if
             
             % calculate the area of the unit cell
diff --git a/helpers/functions/finderb.m b/helpers/functions/finderb.m
index 1bea05a..8c4a958 100644
--- a/helpers/functions/finderb.m
+++ b/helpers/functions/finderb.m
@@ -1,6 +1,6 @@
 %% finderb
 % Binary search algorithm for sorted lists.
-% Searches for the first index _i_ of list where _key_ >= _list(i)_.
+% Searches for the first index _i_ of list where _key_ <= _list(i)_.
 % _key_ can be a scalar or a vector of keys. 
 % _list_ must be a sorted vector.
 % author: André Bojahr