> 
restart: with(Optimization): with(LinearAlgebra): macro(LA= LinearAlgebra): L:=1: r:=2: tau:= 1: interface(rtablesize= 2*r+1): Z:= Matrix( 2*r+1, 2*r+1, [tau, seq(evalf((L/(2*(iz1)*Pi))*sin(2*(iz1)*Pi*tau/L)), iz= 2..r+1), seq(evalf((L/(2*(iz1r)*Pi))*(1cos(2*(iz1r)*Pi*tau/L))), iz= r+2..2*r+1) ], scan= columns, datatype= float[8] ); Dtau00:= < 1 >: Dtau01:= Vector[row](r): Dtau02:= Vector[row](r): Dtau10:= Vector(r): Dtau20:= Vector(r): Dtau1:= LA:DiagonalMatrix([seq(evalf(cos(2*i*Pi*tau/L)), i= 1..r)]): Dtau2:= LA:DiagonalMatrix([seq(evalf(sin(2*i*Pi*tau/L)), i= 1..r)]): Dtau3:= Dtau2: Dtau4:= copy(Dtau1): Dtau:= < < Dtau00  Dtau01  Dtau02 >, < Dtau10  Dtau1  Dtau2 >, < Dtau20  Dtau3  Dtau4 > >; P00:= < L/2 >: P01:= Vector[row](r): P02:= Vector[row](r, j> evalf(L/j/Pi), datatype= float[8]): P10:= Vector(r): P20:= Vector(r, i> evalf(L/2/i/Pi)): P1:= Matrix(r,r): P2:= LA:DiagonalMatrix(P20): P3:= LA:DiagonalMatrix(P20): P4:= Matrix(r,r): P:= < < P00  P01  P02 >, < P10  P1  P2 >, < P20  P3  P4 > >; interface(rtablesize=2*r+1): # optionally J:=Vector([L, L/2 $ 2*r]): # Matrix([[...]]) would also work here E:=DiagonalMatrix(J); X:= Vector[row](2*r+1,symbol=a); U:=Vector[row](2*r+1,symbol=b); X0:= Vector[row](2*r+1,[1]); G:=Vector[row](2*r+1,[1]); C:=simplify(XX0G.ZX.Dtau.P+X.Dtau.ZU.P);
