function A = applyRobinBoundaryConditions(mesh, pGroup, beta, eps )

if isa(beta, 'function_handle')
  x = mesh.P(1,:);
  y = mesh.P(2,:);
  temp = beta(x,y);
  beta = temp;
elseif length(beta) == 1
  beta = beta*ones(1,mesh.NuDoF);
elseif length(beta) ~= mesh.NuDoF
  error('beta does not fit in size or type');
end
if length(eps) == 1
  temp=eps;
  eps = ones(mesh.NuTri,1);
  eps = temp*eps;
  clear 'temp';
elseif length(eps) == mesh.NuDoF
  temp = (eps(tri(:,1)) + eps(tri(:,2)) + eps(tri(:,3)))./3;
  eps = temp;
  clear 'temp';
end

A = sparse(mesh.NuDoF,mesh.NuDoF);

accessRow = mesh.BoundaryEdges(mesh.GBoundaryEdges == pGroup);
firstNode  = mesh.edges(mesh.BoundaryEdges(:),1);
secondNode = mesh.edges(mesh.BoundaryEdges(:),2);

h(mesh.BoundaryEdges(:)) = (mesh.P(1,firstNode)-mesh.P(1,secondNode) ).^2 + ( mesh.P(2,firstNode)-mesh.P(2,secondNode) ).^2;
h = sqrt(h);

for i=accessRow'
  myeps = eps(mesh.edge2triangle(i));
  node1 = mesh.edges(i,1);
  node2 = mesh.edges(i,2);
  A(node1,node1) = A(node1,node1) + h(i)*myeps/12*(3*beta(node1)+beta(node2));
  A(node2,node2) = A(node2,node2) + h(i)*myeps/12*(3*beta(node2)+beta(node1));
  A(node1,node2) = A(node1,node2) + h(i)*myeps/12*(beta(node1)+beta(node2));
  A(node2,node1) = A(node2,node1) + h(i)*myeps/12*(beta(node1)+beta(node2));
end
