diff --git a/classes/heat.m b/classes/heat.m
index 3e1a2bb..0090aa2 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,7 @@ 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));            
             % these are the parameters of the differential equation as they
             % are defined in matlab for the pdesolver
             
@@ -731,11 +731,10 @@ 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;
 
             for k=1:obj.S.numSubSystems
-                parfor i=1:size(tempMap,1)
-                    for n=1:Cells
+                for i=1:size(tempMap,1)
+                    for n=1:obj.S.getNumberOfUnitCells
                         energyMap(i,n,k) = UCmasses(n)*( intHeatCapacity{n,k}(tempMap(i,n,k)) - intHeatCapacity{n,k}(initTemp(n,k)) );
                     end
                 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