A More Efficient Implementation
> x := x0; # initial population
> P := evalm( x ):
> for i from 1998 to 2020 do
> x := evalm( A &* x ); # another year passes
> P := augment( P, x );
> od:
> evalm( transpose( P ) );
>
> T := vector( 24, i -> add( P[j,i], j=1..3 ) );
> PT := stackmatrix( P, T ): # include totals as row 4
> evalm( transpose( PT ) );
>