-Can we solve symbolically matrix equations where the matrices are ungiven by MAPLE?

For a simple example;

Let A.B+2B=C.

where A nxn is matrix and B, C and nx1 matrices. (all matrices are ungiven)

Question: Find symbolically matrix B in terms of C and A?

My try:

We can quickly calculate Matrix B by hand.  We have (A+2I)B=C and so B=(A+2I)^(-1).C   ( where (A+2I) is invertible and I is nxn identity matrix)

Briefly: Can we do for the more complex and difficult equations by MAPLE? 

My question has two steps:

STEP 1:  The multiplication  of is defined as follows

if n<>l, then

.

if n=l and m<=s,

Question 1: I wrote a code for calculating the multiplication  of. Is it right?

The code for Step 1  

restart;

multiply:=proc(n,m,l,s) local g,a: 
a:=unapply(doublefactorial(2*j-1)/factorial(j),j):
g:=unapply((a(m-j)*a(j)*a(s-j)/a(m+s-j))*(2*m+2*s-4*j+1)/(2*m+2*s-2*j+1),j):

if n<>l then 0 else
sqrt((2*m+1)*(2*s+1))*2^(K/2-1).add((g(j)/sqrt(m+s-2*j+1/2))*phi[n, m+s-2*j],j=0..m) 
end if
end proc:
 
n:=2:
l:=2:
m:=1:
s:=1:
multiply(n,m,l,s);

when I compared the results which I got and the results which is given in the book as follows, I think it is right.

Step 2:

We know that the outer product matrix is calculated as follows 

We found the elements of the outer product matrix in Step 1. 

Question 2 : I want to write the elements which are derived in step 1 to the outer product matrix in step 2. In here, the outer product matrix is NxN matrix. N=(M+1).2^(K-1) where K, M are any integers.

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

I have a vector A. I want to re-index vector A.

After re-indexing it, I will use the elements of the vector in new calculations.

For example:

k:=2:
M:=3:
A:=Vector[column]([seq(seq(p*q,q=0..M-1),p=1..2^(k-1))]);
C:=Vector[column]([seq(seq(c(p,q),q=0..M-1),p=1..2^(k-1))]); 
Equate(C,A);
c(1,0)+c(2,1);

c(1,0)+c(2,1)=2.

But the above code doesn' t work.

Suppose that we want to solve second order a PDE by separation of variables method.

Let solution is Summation(X(x).T(t)),

Let X(x):= nPix/l;

And we want to find T(t) by numerical methods. After finding T(t) by numerical methods,

How to differentiate to Summation(X(x).T(t)) with respect to t?   

Briefly; I mean that

How to find derivative of u(t) after finding numerically as follows?

restart:
ode:=diff(u(t),t,t)+diff(u(t),t)-6*u(t)=0:
ICs:=u(0)=1,D(u)(0)=0:
dsolve({ode,ICs},numeric);