function [mat,vector] = DiffusionTerm(D, I, J, S, D_Inv, NumInsideNodes,ConstDiffCoef,Sol)


% Diffusion term's spatial discretization using central differences

if nargin<8    % implicit TS assembles the coeff matrix 'mat' (LHS)
  vector = [];
  
  % neighbor contributions 
  if ConstDiffCoef
    mat = D*sparse(I, J, S.*D_Inv, NumInsideNodes, NumInsideNodes);
  else
    mat = sparse(I, J, D.*S.*D_Inv, NumInsideNodes, NumInsideNodes);
  end
  
  % main diagonal contributions
  mat = sparse(1:NumInsideNodes, 1:NumInsideNodes, sum(mat, 2), NumInsideNodes, NumInsideNodes) - mat;
  
else            % explicit TS assembles the 'vector = mat*sol' (RHS)
  mat = [];
  
  % neighbor contributions
  if ConstDiffCoef
    A = D*sparse(I, J, S.*D_Inv, NumInsideNodes, NumInsideNodes);
  else
    A = sparse(I, J, D.*S.*D_Inv, NumInsideNodes, NumInsideNodes);
  end
  
  % main diagonal contributions
  A = sparse(1:NumInsideNodes, 1:NumInsideNodes, -sum(A, 2), NumInsideNodes, NumInsideNodes) + A;
  vector = A*Sol;
    
end

end  % end of function 'DiffusionTerm'