  function [x] = gaussElim(A,b)
 
  %%% az input ellenőrzése: 
  [n, m] = size(A);
  [nb, mb] = size(b);
  x = nan(n,1);
  if (n ~= m) | (mb ~= 1) | (n ~= nb)
    error("A és/vagy b dimenziója nem megfelelő!")
  endif
  
  %%% Gauss elimináció - "helyben" felülírja az A mátrixot és a b vektort! 
   for j=1:n-1                              % a j-dik ismeretlen kiküszöbölése
       for i=j+1:n                          % a j+1-dik egyenlettől kezdve
          if A(j,j) == 0 
            error("A főelem 0")             % ezt szebben kellene ;-)
          endif
          A(i,j) = A(i,j)/A(j,j);           % az i-dik sor eliminációs konstansa
          for k=j+1:n                       % i-dik sorból levonva a j-dik konstansszorosát
              A(i,k)=A(i,k)-A(i,j)*A(j,k);
          endfor
          % A(i,:) = A(i,:) - A(j,:)*l;     % az előző for ciklus helyett tömörebben így is lehetne
          b(i) = b(i) - A(i,j)*b(j);        % a b vektor i-dik sorából levonva a j-dik konstansszorosát
       endfor
   endfor
  
  % Visszahelyettesítés:
  x(n) = b(n) / A(n,n);               % ???!!! mit tudunk erről?
  
  for i= n-1:-1:1
    sum = 0;
    for j=i+1:n
      sum = sum + A(i,j) * x(j);
    end;
    x(i) = (b(i)- sum) / A(i,i);
  end
  
  %%%% Ehelyett lehetne:
  %U = triu(A)                    % kiveszi A felső háromszög részét
  %x = FelsoHaromszogMegold(U,b)
