function [P,L,U] = pludemo(A);
% PLUDEMO 	Demonstrerar LU-faktorisering som matrismultiplikationer.

[m,n] = size(A);
P = eye(m);
L = eye(m);
U = A;
L1 = eye(m);
% Compute the tolerance
tol = max(m,n)*eps*norm(A,'inf');

% Loop over the entire matrix.
i = 1;
j = 1;
while (i <= m-1) & (j <= n)
   % Find value and index of largest element in the remainder 
   % of column j.
     [y,k] = max(abs(U(i:m,j)))
k = k+i-1;
   %    
   if (y <= tol)
      % The column is negligible, zero it out.
      U(i:m,j) = zeros(m-i+1,1);
      j = j + 1;
   else
y=U(k,j);
clc
disp('Ny omgång påbörjas');
U
disp(['pivotelementets plats är ' '(' int2str(k) ','  int2str(j) ')'])
pause

      % put pivot element in i-th row   
      % Swap i-th and k-th rows.

p=eye(m);
p([i k],:)=p([k i],:);
disp(' ');
disp(['Rad ' int2str(i) ' och rad ' int2str(k) ' byter plats genom matrismultiplikationen U=p*U'])
p
U=p*U

% uppdatera P   P håller reda på vilka radbyten som hittils gjorts
disp('P uppdateras,  P=p*P');
P=p*P
pause

      % Subtract multiples of the pivot row from the lower rows.
l=eye(m);
l(i+1:m,i)=-U(i+1:m,j)/y;
clc
disp('U efter radbyte');
U
disp('Elimination genomförs med hjälp av matrismultiplikationen  U=l*U')
l
U=l*U


% uppdatera L
% När vi är klara har vi fått E*A=U, dvs A=inv(E)*U,  med E=l*p*l*p*.....l*p
% Då är P*A=P*inv(E)*U=p*p*...p*p*inv(l)*p*inv(l)*....p*inv(l)*U=L*U.
disp('L uppdateras,  L=p*L*p*inv(l)')
L=p*L*p*inv(l)
% för att slippa invertera l och multiplicera matriser kan man inse 
% att  L=p*L*p*inv(l) innebär detsamma som följande två rader.
% L([i k],1:i-1)=L([k i],1:i-1);
% L(i+1:m,i)=-l(i+1:m,i)
pause
      i = i + 1;
      j = j + 1;
   end
end

clc
disp('Nu är det klart');
P
L
U
pause
clc
disp('Jämförelse mellan differensen P*A-L*U och tol');
differensen=P*A-L*U
tol