I believe the below is correct.
I=eye(4);
u=zeros(4);
u(1,1:4)=a(1,1:4);
I(1:4,1)=a(1:4,1);
for i=2:4
for j=1:4
u(i,j)=a(i,j)-sum(I(i,1:i-1).*u(1:i-1,j)');
if i<j
I(j,i)=(I(j,i)-sum(I(j,1:i-1).*u(1:i-1,i)'))/u(i,i);
end
end
end
