%  Maximum time step for explicit integration of finite differencing equations 
function [tauMin tau] = dtMaxFDEqs(Domain,Tref,t)
  if iscell(Tref)
    for j=1:length(Tref)
      [tauMin(j) tau{j}]=TimeConst(Domain{j},Tref,t,size(Tref{j}),j);
    end
    tauMin=min(tauMin);
  else
    [tauMin tau]=TimeConst(Domain,Tref,t,size(Tref),0);
  end
end

function [tauMin tau] = TimeConst(Domain,Tref,t,DimT,j)
  [FDEq,nodes]=Domain{1:2};
  Rate=zeros(DimT);
  tau=Inf*ones(DimT);
  Tpert=Tref;
  Dim=sum(DimT>1);
  for i=1:length(FDEq)
    if Dim<=1                             % 1-D
      if j
        Rate(nodes{i})=FDEq{i}(Tref,t,nodes{i},j);
      else
        Rate(nodes{i})=FDEq{i}(Tref,t,nodes{i});
      end
      for m=nodes{i}
        if j
          Tpert{j}(m)=Tref{j}(m)+1;
          Change=FDEq{i}(Tpert,t,m,j)-Rate(m);
        else
          Tpert(m)=Tref(m)+1;
          Change=FDEq{i}(Tpert,t,m)-Rate(m);
        end
        if Change tau(m)=abs(1/Change);
        end
        if j
          Tpert{j}(m)=Tref{j}(m);
        else
          Tpert(m)=Tref(m);
        end
      end
    elseif Dim==2                         % 2-D
      if j
        Rate(nodes{i,1},nodes{i,2})=FDEq{i}(Tref,t,nodes{i,1},nodes{i,2},j);
      else
        Rate(nodes{i,1},nodes{i,2})=FDEq{i}(Tref,t,nodes{i,1},nodes{i,2});
      end
      for m=nodes{i,1}
        for n=nodes{i,2}
          if j
            Tpert{j}(m,n)=Tref{j}(m,n)+1;
            Change=FDEq{i}(Tpert,t,m,n,j)-Rate(m,n);
          else
            Tpert(m,n)=Tref(m,n)+1;
            Change=FDEq{i}(Tpert,t,m,n)-Rate(m,n);
          end
          if Change tau(m,n)=abs(1/Change);
          end
          if j
            Tpert{j}(m,n)=Tref{j}(m,n);
          else
            Tpert(m,n)=Tref(m,n);
          end
        end
      end
    elseif Dim==3                         % 3-D
      if j
        Rate(nodes{i,1},nodes{i,2},nodes{i,3})= ...
                       FDEq{i}(Tref,t,nodes{i,1},nodes{i,2},nodes{i,3},j);
      else
        Rate(nodes{i,1},nodes{i,2},nodes{i,3})= ...
                       FDEq{i}(Tref,t,nodes{i,1},nodes{i,2},nodes{i,3});
      end
      for m=nodes{i,1}
        for n=nodes{i,2}
          for o=nodes{i,3}
            if j
              Tpert{j}(m,n,o)=Tref{j}(m,n,o)+1;
              Change=FDEq{i}(Tpert,t,m,n,o,j)-Rate(m,n,o);
            else
              Tpert(m,n,o)=Tref(m,n,o)+1;
              Change=FDEq{i}(Tpert,t,m,n,o)-Rate(m,n,o);
            end
            if Change tau(m,n,o)=abs(1/Change);
            end
            if j
              Tpert{j}(m,n,o)=Tref{j}(m,n,o);
            else
              Tpert(m,n,o)=Tref(m,n,o);
            end
          end
        end
      end
    end
  end
  tauMin=min(min(min(tau)));
end
