'NULLSPACE - RREF Command in Matlab bugs

Can someone help me debug this Matlab code? Like I have a matrix A, and I need to find A^k by diagonalization method. Below is my original code, but there is a problem with this part:

m = (M - R(i,1)*eye(NumRowsR));
disp(rref(m));
t = null(rref(m));
disp(t); 

This part can't give me the nullspace of matrix t after reduced for some reasons (I recently did some research on the Internet and see some people said about the bug of rref and null). The problem is that it keeps showing me elementary matrix

1 0 0
0 1 0
0 0 1

for the eigenvalues. Does anyone know how to fix it? I will be very much appreciated! Below is my original code, for a better view:

syms x NumRowsM NumColsM NumRowsR NumColsR NumRowst NumColst k numeigen
M = input('Input a square matrix M: ');
k = input('Input the power of matrix M: ');

[NumRowsM, NumColsM]=size(M);  
if (NumRowsM ~= NumColsM)
    disp('Not valid input. Matrix must be a square matrix!')
    
else 
    %Find eigenvalues:
    R = solve(det(M - x*eye(NumRowsM)), x);
    
    [NumRowsR, NumColsR] = size(R);
    if (or(NumRowsR == 0,NumColsR == 0))
        disp('No eigenvalues. The matrix is not diagonalizable')

    else
        numeigen = 0;
        F = zeros(NumRowsR, NumRowsR);
        d = zeros(NumRowsR,1);
    
        for i = 1:NumRowsR
            m = (M - R(i,1)*eye(NumRowsR));
            disp(rref(m));
            t = null(rref(m));
            disp(t);
            [NumRowst, NumColst] = size(t);
        
            if (NumColst == 0)
                if (i == NumRowsR && numeigen > NumRowsR)
                    disp('Matrix not diagonalizable due to invalid eigenvalue');
                    return

                else
                    continue
                end
            else
                numeigen = numeigen + 1;
                if (NumColst == 1)
                    for j = 1:NumRowsR
                        [n, d(j)] = numden(sym(t(j,1)));
                    end
                    for j = 1:NumRowsR
                        F(j,i) = sym(t(j,1)) * lcm(sym(d));
                    end
                else
                    for k = 1:NumColst
                        for j = 1:NumRowsR
                            [n, d(j)] = numden(sym(t(j,k)));
                        end
                        for j = 1:NumRowsR
                                F(j,k) = sym(t(j,k)) * lcm(sym(d));
                        end
                    end
                end
            end
        end
    
    disp(F);
    D = (F\M)*F;
    disp('The power k of the matrix M is: ');
    T = F*(D^k)/(F);
    disp(T)
    end
end


Sources

This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.

Source: Stack Overflow

Solution Source