function [mat_x,vector_x,mat_y,vector_y] = GradientTerm(NatNeighborPairs, PosX, PosY, I, J, S, NumInsideNodes, Sol)

% Check the input arguments to see whether we need to assemble the system
% matrices or the resulting vectors.

if nargin<8    % assembles the coeff matrices 'mat_x' , 'mat_y' (LHS)
  vector_x = []; vector_y = [];
  
  AdvectionField = {1 0};
  [EdgeVelocities_x] = AdvectionVelocity(AdvectionField, NatNeighborPairs, PosX, PosY, I, J, NumInsideNodes, 'CentralDifference');

  AdvectionField = {0 1};
  [EdgeVelocities_y] = AdvectionVelocity(AdvectionField, NatNeighborPairs, PosX, PosY, I, J, NumInsideNodes, 'CentralDifference');  
  
  % neighbor contributions
  mat_x = 1/2*sparse(I, J, EdgeVelocities_x.*S, NumInsideNodes,NumInsideNodes);
  mat_y = 1/2*sparse(I, J, EdgeVelocities_y.*S, NumInsideNodes,NumInsideNodes);
  
  % gradient matrix
  mat_x = sparse(1:NumInsideNodes, 1:NumInsideNodes, sum(mat_x, 2), NumInsideNodes, NumInsideNodes) + mat_x;
  mat_y = sparse(1:NumInsideNodes, 1:NumInsideNodes, sum(mat_y, 2), NumInsideNodes, NumInsideNodes) + mat_y;
  
else            % assembles the 'vector_x = mat_x*sol' (RHS)
  mat_x = []; mat_y = [];
  
  AdvectionField = {@(sol) sol @(sol) 0};
  [EdgeVelocities_x] = AdvectionVelocity(AdvectionField, NatNeighborPairs, PosX, PosY, I, J, NumInsideNodes, 'CentralDifference', Sol);

  AdvectionField = {@(sol) 0 @(sol) sol};
  [EdgeVelocities_y] = AdvectionVelocity(AdvectionField, NatNeighborPairs, PosX, PosY, I, J, NumInsideNodes, 'CentralDifference', Sol);
  
  % neighbor contributions 
  A_x = sparse(I, J, EdgeVelocities_x.*S, NumInsideNodes,NumInsideNodes);
  A_y = sparse(I, J, EdgeVelocities_y.*S, NumInsideNodes,NumInsideNodes);
  
  % gradient vector
  vector_x = sum(A_x, 2);
  vector_y = sum(A_y, 2);
  
end

end  % end of function 'GradientTerm'