> restart; with(linalg):
#  A farmer starts a breeding program where he always fertilizes with a
# plant of genotype dd.  How will the genotype frequencies change over
# time?  The transition matrix is
> A := matrix([[1, 1/2, 0], [0, 1/2, 1], [0,0,0] ]);
# Compute its eigenvalues:
> p := charpoly(A, x);
> lambda := solve(p = 0);
# Next compute the eigenspaces.  We need to solve the equation
# (A-lambda*I)x = 0. 
> R := rref(A - lambda[1]*diag(1,1,1)); 
> evalm(R&*vector([alpha, beta, gamma]));
# So A - lambda1*Id sends the vector (alpha, beta, gamma) to (0,0,0) if
# and only if beta = gamma = 0.  So the eigenspace is just all vectors
# of the form (alpha, 0, 0).  This eigenspace has basis (1,0,0).  We can
# make Maple do the computation for us:
> eigenspace := nullspace(A - lambda[1]*diag(1,1,1));
> v1 := eigenspace[1];
# Repeat for the next eigenvalue, lambda[2] = 1/2:
> R := rref(A - lambda[2]*diag(1,1,1)); 
> evalm(R&*vector([alpha, beta, gamma]));
# The matrix R sends (alpha, beta, gamma) to (0,0,0) if and only if
# alpha + beta = 0 and gamma = 0.  So the eigenspace is all vectors of
# the form (-beta, beta, 0) and it has basis (-1,1,0).
> eigenspace := nullspace(A - lambda[2]*diag(1,1,1));
> v2 := eigenspace[1];
# Repeat for eigenvalue lambda[3] = 0:
> R := rref(A - lambda[3]*diag(1,1,1)); 
> evalm(R&*vector([alpha, beta, gamma]));
# So R sends (alpha, beta, gamma) to (0,0,0) if and only if alpha =
# gamma and beta = -2gamma.  Thus the eigenspace is all vectors of the
# form (gamma, -2gamma, gamma) and it has basis (1,-2,1).
> eigenspace := nullspace(A - lambda[3]*diag(1,1,1));
> v3 := eigenspace[1];
# Form the matrix P whose columns are eigenvectors for A:
> P := augment(v1,v2,v3);
> inverse(P);
# Make a diagonal matrix whose entries are the eigenvalues of A:
> DD := diag(lambda[1], lambda[2], lambda[3]);
> evalm(P&*DD&*inverse(P));
> DDk := diag(lambda[1]^k, lambda[2]^k, lambda[3]^k);
> evalm(P &* DDk &* inverse(P));
# As k goes to infinity this will converge to the matrix
> M := matrix([ [1,1,1], [0,0,0],[0,0,0]]);
# and the equilibrium state will be
> evalm(M &* vector([a_0,b_0,c_0]));
# Example 2:  Suppose the farmer fertilizes each plant with another
# plant of the same genotype.  Then the transition matrix will be
> A := matrix([ [1, 1/4, 0], [0,1/2,0], [0,1/4, 1]]);
> p := charpoly(A, x);
> lambda := solve(p = 0);
# NOTE THAT IN THIS CASE WE HAVE A REPEATED EIGENVALUE.
> R := rref(A - lambda[1]*diag(1,1,1)); 
> evalm(R&*vector([alpha, beta, gamma]));
> eigenspace := nullspace(A - lambda[1]*diag(1,1,1));
> v1 := eigenspace[1];
> v2 := eigenspace[2];
> R := rref(A - lambda[3]*diag(1,1,1)); 
> evalm(R&*vector([alpha, beta, gamma]));
> eigenspace := nullspace(A - lambda[3]*diag(1,1,1));
> v3 := eigenspace[1];
> P := augment(v1,v2,v3); evalm(inverse(P));
> DD := diag(lambda[1], lambda[2], lambda[3]);
> evalm(P&*DD&*inverse(P));
> DDk := diag(lambda[1]^k, lambda[2]^k, lambda[3]^k);
> evalm(P &* DDk &* inverse(P));
# In the limit as k goes to infinity,  we have the matrix
> M := matrix([ [1, 1/2,0], [0,0,0], [0,1/2,1]]);
# and equilibrium state
> evalm(M&*vector([a0,b0,c0]));