Implement a program in Matlab for LU decomposition with pivoting I need to write a program to solve matrix equations Ax=b where A is an nxn matrix, and b is a vector with n entries using LU decomposition. Unfortunately I'm not allowed to use any prewritten codes in Matlab. I am having problems with the first part of my code where i decompose the matrix in to an upper and lower matrix. 
I have written the following code, and cannot work out why it is giving me all zeros on the diagonal in my lower matrix. Any improvements would be greatly appreciated.
function[L R]=LR2(A)

%Decomposition of Matrix AA: A = L R

z=size(A,1);

L=zeros(z,z);

R=zeros(z,z);

for i=1:z

% Finding L

for k=1:i-1

L(i,k)=A(i,k);

for j=1:k-1

L(i,k)= L(i,k)-L(i,j)*R(j,k);

end

L(i,k) = L(i,k)/R(k,k);

end

% Finding R

for k=i:z

R(i,k) = A(i,k);

for j=1:i-1

R(i,k)= R(i,k)-L(i,j)*R(j,k);

end

end

end

R

L

end

I know that i could simply assign all diagonal components of L to be 1, but would like to understand what the problem is with my program! 
I am also wondering how to change this program to include pivoting. I understand I need to say that if a diagonal element is equal to zero something needs to be changed. How would I go about this? 
Thanks in advance!
 A: For backward and forward elimination I used

% Now use a vector y to solve 'Ly=b' 
 for j=1:z

    for k=1:j-1
       b(j)=b(j)-L(j,k)*b(k);
    end;
    b(j) =b(j)/L(j,j); 
 end;

% Now we use this y to solve Rx = y 

x = zeros( z, 1 ); 
for i=z:-1:1    
 x(i) = ( b(i) - R(i, :)*x )/R(i, i); 
end


I put that into your code, and it works ;)
For helping you with pivot strategies it would be helpful to know what strategie you want/have to use as there are various versions. One simple version would be to just swap rows such that the diagonal element $a_{ii} \neq 0$ for all $i = 1, \dots, z$. 
A: 
function[L R x]=LR2(A,b)
% This program will find a solution to
  Ax=b using first giving the
  decomposition of the matrix into L and
  R and then solving.
% Part 1 - Is this matrix square and
  nonsingular? give an error if not
[z y]=size(A);
if (z ~= y )
disp ( 'LR2 error: Matrix must be
  square' );
return;
end;
if det(A)==0
disp('L singular error');

return;

end
% Part 2 : Decomposition of matrix
  into L and R
L=zeros(z,z);
R=zeros(z,z);
for i=1:z
% Finding L
for k=1:i-1
L(i,k)=A(i,k);
for j=1:k-1
L(i,k)= L(i,k)-L(i,j)*R(j,k);
end
L(i,k) = L(i,k)/R(k,k);
end
% Finding R
for k=i:z
R(i,k) = A(i,k);
for j=1:i-1
R(i,k)= R(i,k)-L(i,j)*R(j,k);
end
end
end
for i=1:z
L(i,i)=1;

end
% Program shows R and L
R
L
% Now use a vector y to solve 'Ly=b'
y=zeros(z,1);
y(1)=b(1)/L(1,1);
for i=2:z
y(i)=-L(i,1)*y(1);

for k=2:i-1

    y(i)=y(i)-L(i,k)*y(k);

    y(i)=(b(i)+y(i))/L(i,i);

end;

end;
% Now we use this y to solve Rx = y
x=zeros(z,1);
x(1)=y(1)/R(1,1);
for i=2:z
x(i)=-R(i,1)*x(1);

for k=i:z

    x(i)=x(i)-R(i,k)*x(k);

    x(i)=(y(i)+x(i))/R(i,i);

end;

 end

%print x 
  x 
  end

I have solved the original question. However there are still some problems with my solving of Ly=b and Rx=y. I also need to include pivoting. Any improvements would be greatly appreciated again!
A: After finding the L and U matrix, I tried to find the appropriate B matrix. After then, finding X is a little bit more easier than the traditional way.
% Modifying B
B_new = B;
for i = 2:n
   for k = 1:i-1
       B_new(i,1) = B_new(i,1) - B_new(k,1)*L(i,k); 
   end
end

% Finding X
x = zeros(n, 1);
i = n;
while(i>0)
   x(i) = B_new(i,1);
   k = n;
   while(k > i)
       x(i) = x(i) - x(k)*U(i,k);
       k = k - 1;
   end
   x(i) = x(i) / U(i,i);
   i = i - 1;
end

x

I hope it works for you :)
Greetings from Turkey
NOTE: n equals to z, B equals to b 
