ext_euclid(a,b):=block(
    [x,y,d,x_old,y_old,d_old],
    if b = 0 then return([1,0,a])
    else ([x_old,y_old,d_old]:ext_euclid(b,mod(a,b)),
          [x,y,d]:[y_old,x_old-quotient(a,b)*y_old,d_old],
          return([x,y,d])))$

rref(a):=block([p,q,k],[p,q]:matrix_size(a),a:echelon(a),
    k:min(p,q),
    for i thru min(p,q) do (if a[i,i]=0 then (k:i-1,return())),
    for i:k thru 2 step -1 do (for j from i-1 thru 1 step -1 do a:rowop(a,j,i,a[j,i])),
    a)$

xxx:transpose(matrix([x,y]))$
I:diagmatrix(2,1)$

A:matrix([-2,3],[6,1]);

Eigenvalues

A-r*I;
determinant(%)$
ratsimp(%);
sol:solve(%,r);

Eigenvectors

substitute(sol[1],A-r*I);
rref(%);
row(%,1).xxx;
solve(%,y)$
u:substitute(append(%,[x=1]),xxx);

substitute(sol[2],A-r*I);
rref(%);
row(%,1).xxx;
solve(%,y)$
v:substitute(append(%,[x=1]),xxx);

Diagonalize

S:addcol(u,v);
A.S;
invert(S);
%.A.S;

Fundamental solution for x'(t)=A.x(t)

matrix([substitute(sol[1],exp(r*t)),0],[0,substitute(sol[2],exp(r*t))]);
F:S.%;

Check for linear independence

determinant(F);

General solution

g:F.transpose(matrix([c1,c2]));

Check

ratsimp(diff(g,t)-A.g);

Using initial conditions

invert(S).xxx;
p:F.%;

Check equation and initial condition

ratsimp(diff(p,t)-A.p);
ratsimp(substitute(t=0,p));

Plotdf

tt:A.xxx$
rhs:[tt[1][1],tt[2][1]];
plotdf(rhs, [x,y],
        [tstep,0.01],
        [x,-2,2], [y,-2,2],
        [xfun,"2*x;-x"],
        [direction,forward],
        [nsteps,900])$