function N = normal_derivative_matrix_3D(V,F)
  
  n = size(V,1);
  
  C = cotangent(V,F);
%  C = sym('c', [size(F,1),6]);
  
  [~,on_b] = on_boundary(F);
  [on_b_f,on_b_c] = find(on_b);
    
  
  N=sparse(n, n);
%  N= sym('n',[n, n]);
%  N = 0*N;
  
  reverse_edge = [
  0, 3, 2, 4; 
  3, 0, 1, 5; 
  2, 1, 0, 6; 
  4, 5, 6, 0];

  nb_tets_on_b = size(on_b_f, 1);
  for i=1:nb_tets_on_b
  
    # vertices on boundary subscript
    p_i = mod(on_b_c(i)+1,4)+1;
    q_i = mod(on_b_c(i)+2,4)+1;
    r_i = mod(on_b_c(i)+0,4)+1;
    
    # vertices on boundary indices
    p = F(on_b_f(i), p_i);
    q = F(on_b_f(i), q_i);
    r = F(on_b_f(i), r_i);
    
    # opposite vertex:
    s_i = on_b_c(i);
    s = F(on_b_f(i), s_i);
    
    # collecting weights
    Cpq = C(on_b_f(i), reverse_edge(r_i, s_i));
    Cqr = C(on_b_f(i), reverse_edge(p_i, s_i));
    Crp = C(on_b_f(i), reverse_edge(q_i, s_i));
    
    # setting weights
    N(p, p) = N(p, p) + Cqr;
    N(p, s) = N(p, s) - Cqr;
    
    N(p, q) = N(p, q) + Crp;
    N(p, s) = N(p, s) - Crp;
    
    N(p, r) = N(p, r) + Cpq;
    N(p, s) = N(p, s) - Cpq;
    
    
    N(q, p) = N(q, p) + Cqr;
    N(q, s) = N(q, s) - Cqr;
    
    N(q, q) = N(q, q) + Crp;
    N(q, s) = N(q, s) - Crp;
    
    N(q, r) = N(q, r) + Cpq;
    N(q, s) = N(q, s) - Cpq;
    
    
    N(r, p) = N(r, p) + Cqr;
    N(r, s) = N(r, s) - Cqr;
    
    N(r, q) = N(r, q) + Crp;
    N(r, s) = N(r, s) - Crp;
    
    N(r, r) = N(r, r) + Cpq;
    N(r, s) = N(r, s) - Cpq;
    
  endfor

end