As usual start Maple and load all the linear algebra programs using 
> with(LinearAlgebra);

Create a random 4x4 matrix M:
> randomise(1000*time[real]()); M := RandomMatrix(4, 4, generator = rand(-3 .. 4))
and set two values to unknowns:
> M[2,3]:=x; M[1,4]:=y;
Note that Maple doesn't report the new M unless you ask it to like this:
> M;

Use Maple to calculate the determinant:
> dM:=Determinant(M);
Make a guess at which value of y will guarantee M is non-singular, and check:
> subs(y= number, dM);

Once you get an expression without any x, now look for a pair of values for x 
and y which will make M singular (determinant zero) and substitute them:
> subs({x=number, y=number}, dM);
You should get 0 if you have chosen correctly.


Tell Maple about a new Matrix A:
> A:=<<9, 6, -4>|<-9, -6, 6>|<-3, -3, 5 >>;

Use Maple to check that the matrix is non-singular:
> Determinant(A);
and then form the matrix from which we will get the eigenvalues:
> A1:=A-lambda*IdentityMatrix(3);

Use a determinant row operation then use the maple command Transpose 
on that matrix and do another row operation on the transposed matrix to 
create two zeros. Now use SubMatrix to extract a 2x2 matrix cofactor S 
corresponding to the non-zero element and evaluate
> Determinant(S);
The roots of this polynomial are the other eigenvalues, the root of the 
non-zero element is the other eigenvalue.

Repeat with a different pair of rows to identify a different submatrix determinant
and check that the roots are the same as previously.

We can find the polynomial whose roots are the eigenvalues using:
> pA:=Determinant(A1);
which is the so-called characteristic polynomial of A.
Note that the coefficient of pA without any lambda in is the determinant of A.

To manipulate polynomials you can use factor, expand, sort, simplify, etc. 
> factor(pA);
Verify that all of the different linear polynomials in lambda are included in the factorisation.

You can also plot pA to see the locations of the eigenvalues using
> plot(pA,lambda=1..4);
The values 1 and 4 are the plotting range.

We can also get the system of equations to solve for the eigenvalue equal to 2 using this matrix:
> Aj:=<subs(lambda=2,A1)|<0,0,0>>;
Pivot this matrix to create a row of zeros and get a solution.

Compare what you just got with this vector solution:
> v1:=LinearSolve(subs(lambda=2,A1),<0,0,0>,free=s);
What value of s gives your solution?

Verify that v1 when multiplied by A on the left gives the expected result (for any s).
> Multiply(A,v1);
> 2*v1;

Find two different eigenvectors of the eigenvalue lambda=3 similarly and check them. 
3 is an eigenvalue of multiplicity 2 since it has two non-parallel vector solutions.

What are the homogeneous solutions if you try to substitute a non eigenvalue 
into A1 instead of 2 or 3?


Try to repeat the row operations procedure for 
> B:=<<32, -36, -18>|<27, -31, -18>|<-6, 6, 2 >>; 
> B1:=B-lambda*IdentityMatrix(3);

Some row operations give two easy zeros but some don't. Find a pair which work.

Now try with this matrix
> C:= <<-6, 2, -3>|<-13, 15, -15>|<-2, 8, -7>>; 
> C1:=C-lambda*IdentityMatrix(3);
but note that each way you try to do operations on C1 no easy cancellation occurs.

Use
> Eigenvalues(C);
to note what kind of factors you should get and try to create one. If after 
5 minutes if you have failed to create two zeros in a row, just evaluate 
> fC:=factor(Determinant(C1));
to reassure yourself the answers Maple gave are indeed correct.

Find one of the three eigenvectors in C using RowOperation and/or Pivot and check for 
it or a multiple of it in the columns of P in

> (v,P):=Eigenvectors(C);

Extract another eigenvector from P and check it has the expected eigenvalue.

Make a diagonal matrix of eigenvalues using

> D1:=DiagonalMatrix(v);

Check that P and D1 multiplied is the same as C and P multiplied in that order.

Create a matrix which has random eigenvalues from w as follows:
> randomise(time[real]()):w:=RandomVector(3,generator=rand(1..9));
> Dw:=DiagonalMatrix(w);
> Q:=RandomMatrix(3,3,generator=rand(-3..4));

We choose Q as a random matrix for the eigenvectors and (if Q is 
non-singular) use diagonalisation to make J from Dw as follows:
> J:=Multiply(Multiply(Q,Dw),MatrixInverse(Q));

Check that the eigenvalues and eigenvectors of J are as you expect.