In this paper, based on the Gauss-Jordan elimination to solve linear systems, the algorithms for finding all the magic matrices of size 4 by 4 are developed and discussed, and the efficiency of the algorithms is tested after complementing the algorithms in MATLAB. The total magic matrices of 4 by 4 can be found with the program; for completeness and explaining the algorithms, the case of 3 by 3 is also included.