Gaussian Reduction with Partial Pivoting in Octave / Matlab

Question

I am trying to get an algorithm from an old math book I have to work properly. The book is Cullen, Introduction to Numerical Linear Algebra. Algorithm 2.3.

function B = gaussian2 ( A )
  B = A;
  %INPUT: A(m,n)
  [m,n] = size(B);
  %m1 = min{n, m - 1}
  m1 = min(n,m - 1);
  %For i = 1,2,...,m1
  for i = 1:m1
    % Choose j so that |a_ji| >= |a_ki| for k = i,...,m.
    for k = i:m
      [a j] = max(B(:,i))
    end
    %Interchange rows i and j
    B([i j],:) = B([j i],:)
    %If a_ii = 0, increment i
    if ( B(i,i) <= eps )
      i = i + 1;
    end
    %For k = i + 1,...,m
    for k = i+1:m
      %T = -a_ki/a_ii
      T = -B(k,i)/B(i,i);
      %a_ki = 0
      B(k,i) = 0;
      %For j = i + 1,...,n
      for j = i+1:n
        %a_kj = T*a_ij + a_kj
        B(k,j) = T*B(i,j) + B(k,j);
      end
    end
  end
end
% works for matrix B from book but not matrix A
% A = [ 5 7 6 5;7 10 8 7; 6 8 10 9; 5 7 9 10]
% B = [ 4 -1 -1 0 1; -1 4 -1 -1 2; -1 -1 4 -1 3; 0 -1 -1 4 4]

The comments are the algorithm steps and I believe that my error is in the step 4.

% Choose j so that |a_ji| >= |a_ki| for k = i,...,m.
for k = i:m
  [a j] = max(B(:,i))
end

If someone has an idea on how to implement this properly, or a fix suggestion, please let me know.

The answer matricies are

A = [7 10 8 7; 0 -4/7 22/7 3; 0 0 5/2 17/4; 0 0 0 1/10] 
B = [4 -1 -1 0 1; 0 15/4 -5/4 -1 9/4; 0 0 10/3 -4/3 4; 0 0 0 16/5 31/5]

The current algorithm works with B but fails with A.

function B = gaussian2 ( A )
  %INPUT: A(m,n)
  %B = A
  %m1 = min{n, m - 1}
  %For i = 1,2,...,m1
  %  Choose j so that |B_ji| >= |B_ki| for k = i,...,m.
  %  Interchange rows i and j
  %  If B_ii = 0, increment i
  %  For k = i + 1,...,m
  %    B_ki = 0
  %    For j = i + 1,...,n
  %      B_jp = T*B_ij + B_kj
  B = A;
  [m, n] = size(B);
  m1 = min(n,m - 1);
  for i = 1:m1
    j = max(abs(B(i:m,i)));
    for k = i:m
      if(abs(B(k,i))==j)
        j = k;
        break;
      end
    end
    B([i j],:) = B([j i],:);
    if ( abs(B(i,i)) <= eps )
      i = i + 1;
    end
    for k = i+1:m
      T = -B(k,i)/B(i,i);
      B(k,i) = 0;
      for j = i+1:n
        B(k,j) = T*B(i,j) + B(k,j);
      end
    end
  end
end

This algorithm works with both A and B matricies and gives the proper answers per the book.