Algorithm 2: using orthogonality
First I looked into the matrix pattern in 2x2 and 4x4. I saw the rows are orthogonal. I thought, ``Aha, because the determinant is volume and when a simplex has the maximal volume when the edge vector length is fixed? Orthogonal vectors!'' This is quite intuitive for me.
Therefore, I implemented a method that looked up the orthogonal vectors. This is program 2.
Program 2
function algo_02(mat_rank)
% Introduction to linear algebra Chapter 5. problem 33
% Algorithm 2: generate Hadamard matrix (each row is orthogonal), but
% this only can gives me 1,2,4k matrices
% @author Hitoshi
if nargin ~= 1;
error('Usage: algo_02(mat_rank).')
end
% possible element set A_i = {1, 1}
SetA = [1 1];
cand_mat = zeros(mat_rank, mat_rank);
cand_mat(1,:) = ones(1, mat_rank);
cand_row = zeros(1, mat_rank);
global MAXDET
global MAXDET_MAT
MAXDET = 0;
MAXDET_MAT = zeros(1, mat_rank * mat_rank);
cur_row_index = 2;
loopdepth = 1;
gen_comb_set(SetA, cand_mat, cand_row, mat_rank, loopdepth, cur_row_index);
MAXDET_MAT
fprintf(1, 'max detderminant = %d.\n', MAXDET);
end
%%
% Looking for the orthogonal rows and compute the determinant.
% \param SetA element candidate set
% \param cand_mat current candidate matrix
% \param cand_row current candidate row
% \param mat_rank rank of matrix (not exactly the rank, size of n)
% \param loopdepth parameter to simulate forloop depth by recursion.
% \param cur_row current row index to look for
function gen_comb_set(SetA, cand_mat, cand_row, mat_rank, loopdepth, cur_row)
global MAXDET;
global MAXDET_MAT;
num_set = mat_rank;
num_cand = size(SetA);
szSetA = size(SetA);
% This should be assert(sum(szSetA == [1 2]) == 2)
if sum(szSetA == [1 2]) ~= 2
error('Not assumed set candidate matrix (should be 1x2)')
end
if cur_row > mat_rank;
% cand_mat;
det_a = det(cand_mat);
if det_a > MAXDET
MAXDET = det_a;
MAXDET_MAT = cand_mat;
end
elseif loopdepth > num_set
if check_orthogonal_row(cand_mat, cand_row, cur_row) == 1
cand_mat(cur_row, :) = cand_row;
cand_row = zeros(1, mat_rank);
cur_row = cur_row + 1;
gen_comb_set(SetA, cand_mat, cand_row, mat_rank, 1, cur_row);
end
else
% raw is not yet ready, generate it.
for j = 1:szSetA(2)
cand_row(loopdepth) = SetA(j);
gen_comb_set(SetA, cand_mat, cand_row, mat_rank, loopdepth + ...
1, cur_row);
end
end
end
%%
% check the rows are orthogonal with rows < cur_row
% \param cand_mat current candidate matrix
% \param cand_row current candidate row
% \param cur_row current row index to look for the orthogonal
function ret_is_all_orthogonal = check_orthogonal_row(cand_mat, cand_row, cur_row)
is_all_orthogonal = 1;
for i = 1:(cur_row  1)
if dot(cand_mat(i,:), cand_row) ~= 0
is_all_orthogonal = 0;
break
end
end
ret_is_all_orthogonal = is_all_orthogonal;
end
But, this program gives me the max determinant value is zero when 3x3, 5x5, and 6x6 matrix. This is strange. For instance, I can easily find a non zero determinant matrix, for instance, [1 1 1; 1 1 1 ; 1 1 1] for 3x3. The determinant is 4. Also I realized I can not make a orthogonal rows in 3x3 case as the following.
When I think about the geometry, it is also easy to see it is not possible. Figure 1 shows we can not generates orthogonal vectors in 3D case when the coordinates value are only allowed {1, 1}.

Figure 1: 3D, can not make orthogonal vectors by using {1,1} coordinates 
Figure 2 shows that this method works in 2D case. There are cases that we could make orthogonal vectors even the coordinate values are limited. Figure 2 also shows that the volume (= area, in 2D) that represents the determinant. This is (√2)^2 = 2, and the max determinant of 2x2 matrix is also 2.

Figure 2: Orthogonal vectors by using {1,1} coordinates in 2D case. 
This method can be applied to only 1,2,4n (n >=1) cases. At this point, I found these kind of matrices are called Hadamard's matrix. This problem is called Hadamard's Maximum Determinant Problem. On the Web, there is even the number (max determinant value) for 6x6 case. I am surprised that a lot of cases are known. In the case of 1,2,4n, there is a construction method to generate a Hadamard matrix. The number 1,2,4n is called Hadamard number.
Moreover, matlab/octave has function hadamard(), this generates a Hadamard matrix.
But, I didn't know how to compute the max determinant value of nonHadamard number matrix. According to
http://mathworld.wolfram.com/HadamardsMaximumDeterminantProblem.html, the max determinant value sequence of Hadamard matrix is known in 1962. There should be a clever method.