Question: Writing set of differential equations as a matrix differential equation form

Suppose that we have a set of differential equations as follows

and we want to write the above equation as a matrix form  ( for convenience, Let upper limit of summations be 2, not infinity )

QUESTION: How to find the coefficients of Matrix Dif. equations (matrices M, C, K, f, D) by MAPLE? MY TRY: question.mw

We can find the matrices by hand as follows: (But how to derive the matrices directly from the equation in the first picture by Maple)

restart:with(LinearAlgebra):with(plots):interface(rtablesize=36):  
with(Physics): #for KroneckerDelta
   
N:=2:
p:=2:
A:=Matrix(N,1):
B:=Matrix(N,1):
for n from 1 to N do 
A(n,1):=diff(phi[n](t),t,t)+2*gamma*add(sin(n*Pi*v0*tau)*sin(j*v0*Pi*tau)*diff(phi[j](tau),tau,tau),j=1..N)+4*gamma*Pi*v0*add(j*sin(n*Pi*v0*tau)*cos(j*v0*Pi*tau)*diff(phi[j](tau),tau),j=1..N)+n^4*Pi^2*phi[n]-2*gamma*add( (j*Pi*v0^2)*sin(n*Pi*v0*tau)*sin(j*v0*Pi*tau)* phi[j](tau),j=1..N):
end do: 
B:=Matrix(N,1):
for n from 1 to N do 
B(n,1):=gamma*Pi^3/2*sin(n*Pi*v0*tau)+add( U[i](tau)*sin(n*Pi*eta[i]),i=1..p):
end do:
# in order to make ensure, let's look backwards

r:=2:

M:=Matrix(r,r):
C:=Matrix(r,r):
K:=Matrix(r,r):
for n from 1 to r do
for j from 1 to r do  
M(n,j):=KroneckerDelta[n,j]+2*gamma*sin(n*Pi*v0*tau)*sin(j*Pi*v0*tau);
C(n,j):=4*gamma*Pi*v0*j*sin(n*Pi*v0*tau)*cos(j*Pi*v0*tau):
K(n,j):=n^4*Pi^2*KroneckerDelta[n,j]-2*gamma*(j*Pi*v0)^2*sin(n*Pi*v0*tau)*sin(j*Pi*v0*tau):
end do:
end do:
f:=Matrix(r,1):
X:=Matrix(r,1):
U:=Matrix(r,1):
for n from 1 to r do
f(n,1):=gamma*Pi^3/2*sin(n*Pi*v0*tau);
X(n,1):=phi[n](tau):
U(n,1):=u[n](tau):
end do:
DD:=Matrix(r,p):
for n from 1 to r do
for i from 1 to p do 
DD(n,i):=sin(n*Pi*eta[i]);
end do:
end do:  
 M*diff(X,tau,tau)+ C*diff(X,tau)+ K*X=f+DD.U;

quest.mw