The short answer is that it's possible to compute an orthonormal basis of eigenvectors for an orthogonal matrix, but that MATLAB doesn't check for orthogonal matrices in EIG and so just provides an answer for a generic nonsymmetric matrix. For general nonsymmetric matrices, there typically isn't an orthogonal basis of eigenvectors.
However, if you know that your matrix is orthogonal, you can compute an orthogonal eigenbasis as follows:
% Construct an orthogonal matrix with a repeated eigenvalue 1 (if they are
% all different, eig will already return an orthonormal basis of
% eigenvectors)
rng default;
Q = blkdiag(eye(5), orth(randn(5)));
V = orth(randn(10));
Q = V*Q*V';
eig(Q)
% Compute eigenvalues using eig:
[U, D] = eig(Q);
res = norm(Q*U - U*D) % residual shows these are correct eigenvalues and eigenvectors
orthogonality = norm(U'*U - eye(10)) % but the eigenvectors aren't orthogonal
% Instead, compute the Schur decomposition:
[U, S] = schur(Q, 'complex');
S
% since S is a triangular matrix which we know to also be orthogonal up to round-off,
% it must be diagonal up to round-off. Set D to be a diagonal matrix
% containing the same diagonal as S.
D = diag(diag(S));
res = norm(Q*U - U*D) % residual shows these are correct eigenvalues and eigenvectors
orthogonality = norm(U'*U - eye(10)) % and the eigenvectors are also orthogonal
We're unlikely to make EIG try to detect if an input matrix is orthogonal and do this kind of operation, because we can only determine if a matrix is orthogonal up to a tolerance, and would have to make a choice of a cut-off at which we don't count a matrix as orthogonal anymore. These kinds of hidden behaviors have usually been more trouble than they are worth: imagine if you were passing a bunch of orthogonal matrices to EIG in some function, and relying on the eigenvectors returned being orthogonal. If one of them was slightly further from being orthogonal, the eigenvectors could suddenly be returned as not orthogonal at all for this matrix.
