function H = Heffter_systems_n_4(n,k,t)



  v = 2*n*k+t;      % Dimension of the group
  Sym = 1:v/2;      % Symbols of an half set
  % Minus the subgroup of size t
  aux = find(mod(Sym,v/t)==0);
  Sym(aux) = [];

  Triples = zero_sum_triples(Sym,v); % Zero-sum triples in Sym

  % Without loss of generality, we assume that the first row and the first
  % column contain the element 1 (or -1)
  aux1 = Triples(:,1);
  pos = find(aux1==1);
  A1 = Triples(pos,:);
  pos = find(aux1~=1);
  A2 = Triples(pos,:);
  A3 = A2;
  A4 = A2;

  % When two zero sum triples of A2 are compatibles?
  A_int = zeros(size(A2));
  for i = 1:size(A2,1)
    for j = i+1:size(A2,1)
      if length(intersect(A2(i,:), A2(j,:)))==0
        A_int(i,j) = 1;
      endif
    endfor
  endfor
  A_int = A_int+A_int';



  trovati = 0;

  noniso = 0;
  database = [];
  for i1 = 1:size(A1,1)
    b1 = A1(i1,:);                  % First subset of Heffter system
    for i2 = 1:size(A2,1)
      b2 = A2(i2,:);                % Second subset of Heffter system
      G = [b1, b2];
      index2 = find(A_int(i2,:));   % search the third subset between
      % the compatible ones
      if length(unique(G)) ==length(G)
        for i3 = 1:size(A3(index2,:),1)
          b3 = A3(index2(i3),:);    % Third subset of Heffter system
          G = [b1, b2,b3];
          index3 = intersect(index2,find(A_int(index2(i3),:)));
          if length(unique(G)) ==length(G)
            for i4 = 1:size(A4(index3,:),1)
              b4 = A4(index3(i4),:);
              G = [b1, b2,b3,b4];
              index4 = intersect(index3,find(A_int(index3(i4),:)));
              if length(unique(G)) ==length(G)
                G = [b1; b2;b3;b4];
                trovati = trovati+1;
                new_Heffter_System = true;
                if noniso>0      % Is this a new Heffter system, or is it an old one?
                  for ijk = 1:length(database)
                    Old_Heffter = database{ijk};
                    % order the rows based on their smallest element;
                    for row_index = 1:size(G,1)
                      G(row_index,:) = sort(G(row_index,:));
                      Old_Heffter(row_index,:) = sort(Old_Heffter(row_index,:));
                    endfor
                    [~,idG] = sort(G(:,1));
                    [~,idOld] = sort(Old_Heffter(:,1));
                    Old_Heffter = Old_Heffter(idOld,:);
                    G = G(idG,:);

                    % Check if every row contains the same elements:
                    if norm(G - Old_Heffter)==0
                      new_Heffter_System = false;
                    endif
##                    row_same_el = true;
##                    for row_index = 1:size(G,1)
##                      if max(abs(sort(G(row_index,:)) - sort(Old_Heffter(row_index,:))))~=0
##                        row_same_el = false ;
##                      endif
##                    endfor
                  endfor
                  if new_Heffter_System %(row_same_el == false)
                    noniso = noniso+1;
                    database{noniso} = G;
                  endif
                else
                  noniso = noniso+1;
                  database{noniso} = G;
                endif
              endif
            endfor
          endif
        endfor
      endif
    endfor
  endfor

  disp('Heffter systems for these parameters:')
  database % These are all the Heffter systems on the set of
  % symbols Sym


  % Can we find two Heffter Systems inside of database that are
  % orthogonal?

  mm = length(database);
  num_orthogonal = 0;
  H = [];
  for i = 1:mm
    A1 = database{i};      % Fix one Heffter system
    for j = i+1:mm
      A2 = database{j};     % And the second one
      orthogonal = true;
      for k1 = 1:size(A1,1)  % And now check if they are orthogonal
        for k2 = 1:size(A2,1)
          if length(intersect(A1(k1,:),A2(k2,:)))>1
            orthogonal = false;
          endif
        endfor
      endfor
      if orthogonal % If they are orthogonal, construct the array
        Heff_array = isHeffter(A1,A2);
        % Before adding the array to the database, check if it is already there
        isnew = true;
        Heff_array1 = Heff_array(:,1);
        [~,idx1] = sort(Heff_array1);
        Heff_array = Heff_array(idx1,:);

        for kk = 1:length(H)
          C = H{kk};
          C1= C(:,1); % As before, check up to a permutation of the rows
          if max(abs(sort(Heff_array1) - sort(C1))) == 0
            [~,idx2] = sort(C1);
            C = C(idx2,:);
          endif

          if norm(C - Heff_array)==0
            isnew = false;
          endif
        endfor
        if isnew
          num_orthogonal = num_orthogonal + 1;
          H{num_orthogonal} = Heff_array;
        endif
      endif
    endfor
  endfor


