restart;
  A:= Matrix(4, 8, [ [-1, 1-t, 0, 0, 0, 0, 0, 0],
                     [0, 0, 1, 1+t, 0, 0, 0, 0],
                     [0, 0, 0, 0, -1, 1-t, 0, 0],
                     [0, 0, 0, 0, 0, 0, 1, 1+t]
                   ]
             );
  B:= Matrix(4, 8, [ [-1, 1-t1, 0, 0, 0, 0, 0, 0],
                     [0, 0, 1, 1+t1, 0, 0, 0, 0],
                     [0, 0, 0, 0, -1, 1-t2, 0, 0],
                     [0, 0, 0, 0, 0, 0, 1, 1+t2]
                   ]
            );
#
# Define function which generates target matrix from
# supplied matrix
#
  genMat:= M-> Matrix( op(1, A),
                       (i,j)-> subs( t=cat(t, iquo(i-1, 2)+1),
                                     A[i,j]
                                   )
                     );
  C:=genMat(A);
#
# Check B=C
#
  LinearAlgebra:-Equal(B, C);

  restart;
#
# Define function whihc builds required matrix
#
  doMat:= (n::integer)-> Matrix(n+1, n+1, (i,j)->`if`(j=i+1, i,0)):
#
# Examples of use
#
  M:=doMat(4);
  M:=doMat(8);

  restart:
  x:= 12: y:= 46: s:= 0:
  L:= NULL:
  while 0 < y do
        if   type(y, odd)
        then s:= s + x:
             y:= y - 1:
        else x:= 2*x:
             y:= y/2:
        end if;
        L:= L,x,y,s;
  end do:
  L:= [L];

restart;
#Define a table;
T := table(symmetric, [(1,1) = 1,(2,2)=1,(3,3)=1,(4,4)=1,(1,2)=1,(1,3)=1,(1,4)=1,(2,3)=1,(2,4)=1,(3,4)=1]);

type(T,'table');
type(T,'name');
type(T,'symbol');
type(eval(T),'table');

type(op(1,T), 'table');
type(op(1,T), 'name');
type(op(1,T), 'symbol');

  restart;
  f:= x->x^2;
#
# Unit-spaced values
#
  f~([$1..10])[];
#
# For non-unit-spaced values, use
# elementwise operators
#
  f~([$1..10]/~3)[];
  f~([$1..10]*~1.75)[];
#
# Or for some random list of values
#
  L:=[1, y, 2.5, x^2, 3.7]:
  f~(L)[];

  restart;
#
# Define assorted parameters
#
  Hlist:= [0, 5, 6, 10]:
  params:= {a = 7, alpha = 2}:
#
# Defne the ODE system
#
  odeSys:= { diff(f(x), x$3) + 2*GB*a*alpha*f(x) = 0,
             f(0) = 1,
             D(f)(0) = 0,
             f(1) = 0
           }:

#
# Initialize sequence for output plots and output
# data
#
  plts:= NULL:
  dat:= NULL:
#
# Loop through values in Hlist, generating  solutions,
# plots and data
#
  colors:=[red, blue, green, black]:
  for k to numelems(Hlist) do;
      Sol_f[H]:= dsolve
                  ( eval
                    ( odeSys,
                      params union {GB = Hlist[k]}
                    ),
                    numeric
                  );
    #
    # Generate plot of solution, and add to 'plts' sequence
    #
      plts:= plts,
             plots:-odeplot( Sol_f[H],
                             [x, f(x)],
                             x = 0 .. 1,
                             color=colors[k]
                           );
    #
    # Extract plot data as a matrix and add matrix to 'dat' sequence
    #
      dat:= dat,
            plottools:-getdata([plts][-1])[3];
    #
    # Export this data to labelled sheets of an Excel file
    #
    # NB - OP will need to change filename (highlighted)to something
    # appropriate for his/her machine
    #
      ExcelTools:-Export( [dat][-1],
                          "C:/Users/TomLeslie/Desktop/testDat.xlsx",
                          cat("GB=",Hlist[k]),
                          "B$2"
                          );
  end do:
#
# Plot the solutions (just for fun)
#
  plots:-display([plts] );

  restart;
#
# Define a function which iterates a given
# function 'f', 'n' times
#
  gF:= (f,n)-> unapply( (f@@n)(x), x);
#
# Compare 'gF' with builtin function for
# numeric arguemnts
#
  sin(sin(sin(0.25)));
  gF(sin, 3)(0.25);
#
# Compare 'gF' with builtin function for
# symbolic arguemnts
#
  sin(sin(sin(z)));
  gF(sin, 3)(z);
#
# Compare 'gF' with user-defined function
# for numeric arguemnts
#
  f1:= x-> sin(x)+cos(x)+x^2;
  f1(f1(f1(0.25)));
  gF(f1, 3)(0.25);
#
# Compare 'gF' with user-defined function
# for symbolic arguemnts
#
  f1(f1(f1(z)));
  gF(f1, 3)(z);