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 ) );

>