Update matrixfactorizations.mac

Fix "diagonalize". It was not correctly orthogonalising the eigenspaces if the given matrix was symmetric as described.

Author: LukeLongworth

stack/maxima/contrib/matrixfactorizations.mac

@@ -69,10 +69,10 @@ get_Jordan_form(M):= block([jordan_info,J,P],
  * @param[matrix] M a diagonalizable matrix
  * @return[list] A list of two matrices: [P, D] such that D is diagonal and M = P . D . P^^-1. Returns empty list if M is not diagonalizable
  */
-diagonalize(M):= block([P,D],
+ddiagonalize(M):= block([P,D],
   if not(squarep(M)) then return([]),
   [P, D]: get_Jordan_form(M),
-  if ev(zeromatrixp(M - transpose(M)),simp) then P:ev(transpose(apply('matrix, map(lambda([ex], ex/sqrt(ex.ex)), args(transpose(P))))), simp),
+  if ev(zeromatrixp(M - transpose(M)),simp) then P:ev(transpose(apply('matrix, map(lambda([ex], ex/sqrt(ex.ex)), gramschmidt(args(transpose(P)))))), simp),
   if diagp(D) then return([P,D]) else return([])
 );