As usual start Maple and load all the linear algebra progams.

Create a matrix which has the eigenvalues from w:
> randomise():w:=RandomVector(3,generator=rand(1..9));
> Dw:=DiagonalMatrix(w);
> Q:=RandomMatrix(3,3,generator=rand(-3..3));
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.

Maple can create the kth power of a matrix automatically:
> MatrixPower(Dw,k);
Check that when you use k=2 you get the same as
> Multiply(Dw,Dw);

Create a matrix B and its eigenvalues and eigenvectors:
> B:=Matrix([[6, -2, 1], [7, -3, 1], [-4, 2, 1]]);

See what Maple produces for 
> MBk:=MatrixPower(B,k);
and check that the powers of the eigenvalues of B are in there.

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

Multiply P and a diagonal matrix of the powers of the eigenvalues, and 
the inverse of P to get Bk, the kth power of B.

> Bk:=Multiply(Multiply(P,MatrixPower(DiagonalMatrix(v),k)),MatrixInverse(P));

This should be something like MBk, but might be a bit different in some ways. You can see they are the same using 
> simplify(Bk-MBk);

Check that Bk and MBk are the same for small values of k:

> subs(k=2,Bk); subs(k=2,MBk); Multiply(B,B);

Check also for k=-1 and 0.


The matrix underlying the recurrence c(n+1)= 6*c(n) + 7*c(n-1) - 60*c(n-2) is
> C:=<<6,1,0>|<7,0,1>|<-60,0,0>>;

Diagonalise C utilising
> (v,P):=Eigenvectors(C);
and notice that its eigenvectors have the special "powers
of eigenvalues" formula they are supposed to have.

Form Ck as the kth power of C and now suppose that the first three values in the 
sequence are c(0)=8, c(1)=-85 and c(2)=27. 
Multiply Ck by the initial values vector, remembering to use the correct order.

Calculate c(3) and c(4) using the recurrence then use the power matrix to
get c(k) and check with your values with k=3 and 4.
At which point is c(k) below zero for the last time?

Find initial values for the sequence which would give a solution with only 
two different eigenvalue powers in the final solution for c(k).


You can plot points using this command:
> plot([[-1,5],[2,4],[3,-1]],style=point, symbolsize=30,colour=blue);
The matrix underlying this system of equations is
> E2:=<<1,4,9>|<-1,2,3>|<1,1,1>|<5,4,-1>>;

Find the solutions to this system of equations to get the 
coefficients of a quadratic polynomial which you can then plot using:
> plot(your polynomial here,x=-1..3);
Note that the curve passes through all three points.

Repeat with polynomials of degree 3 through these points:
    [-1,3], [-2,6], [4,-12]
the system of equations in this cases is:
> E3:=<<-1,-8,64>|<1,4,16>|<-1,-2,4>|<1,1,1>|<3,6,-12>>;
Solve for a*x^3+b*x^2 + c*x +d using
> LinearSolve(E3,method='none',free=t);
to see that a=-t, b=t, c=10*t -3, d=8t.

You can plot the family of polynomials of degree 3 like this:

> with(plots);
> p3j:=(-t*x^3 + t*x^2+ (10*t-3)*x + 8*t);
> animate( plot, [p3j,x=-2..4], t=-10..10, trace=5, frames=50 );

To animate the plot you click on the image and then press the 
"play" button that is now in the toolbar at the top. Note how every 
cubic curve passes through all 3 specified points.

Pick 5 random points and create the system of equations for 
a straight line through them:
> xv:=RandomVector(5,generator=-4..4);
> F:=< xv | <1,1,1,1,1>>;
> yv:=RandomVector(5,generator=0..7);

Use this to plot your points:
> plot((xv,yv),style=point, symbolsize=30,colour=red);

And 
> FTF:=Multiply(Transpose(F),F);
> FTy:=Multiply(Transpose(F),yv);
Plot the solution to FTF times v = FTy, which is the best fit line,
by using RREF, inverses or by doing 
> LinearSolve(FTF,FTy); 
and repeat with the same xv and yv to plot the best fit quadratic.