I have an ODE. I solved it numerically by Runge Kutta method.

The solution is X(t).

Now, I want to multiply X(t) by sin(y) and to plot in 3D.

How can we do it? 

question.mw

What is the simplest method for adding a color bar for plottings? (a scale that shows the correspondence 
between numbers and colors)

and could you interpret the graphics and the color bar?  I need your valuable comments.

I use Maple 2018.

For example;  

question.mw

I mean

How can we create the following upper triangular matrix?

where k and M are any integers,
 F and L are MxM Matrices  as follows

restart;
with(LinearAlgebra):
k:=2:
M:=3:
F:=Matrix(M,M):
for i from 1 to M do 
  for j from 1 to M do
 F[1,1] := 2;
    F[i,j]:=0
  end do
end do:
F;
L:=Matrix(M,M): 
L:=LinearAlgebra[BandMatrix]([
      [seq(-sqrt(2*i-1)/(2*i-1)*sqrt(2*i-3),i=2..M)], [1,seq(0*i,i=1..M-1)],[seq((sqrt(2*i-3))/((2*i-3)*sqrt(2*i-1)),i=2..M)]]);

I have a PDE.

I am trying to code for solving the pde by numerically and analytically.

But the code has an error.

Could you help me?

Best regards.
THE CODE: pde.mw

How to solve the following second order ODE system? (I want to solve analytically and numerically and to compare the results) Thanks. Best regards.

where

Initial conditions:

Other symbols in the equation (v_0,a_0 etc.) are constant.

MY CODE TRY:  Code.mw

I tried to write a code for matrices in the question. ( I used x instead of gamma in the code)

restart:
n:=2:
  M:= tau -> Matrix
      ( n,
        n,
        shape=identity
      )
      +Matrix
      ( n,
        n,             (i,j)->2*x*sin(i*Pi*v_0*tau)*sin(j*Pi*v_0*tau)
            ):
  C:=tau -> Matrix
              ( n,
                n,
                (i,j)->4*(i*Pi*v_0*tau)*sin(i*Pi*v_0*tau)*cos(j*Pi*v_0*tau)
              ):
  K:=tau -> n^4*Pi^2* Matrix
      ( n,
        n,
        shape=identity
      )-Matrix
      ( n,
        n,
        (i,j)-> 2*mu*(i*Pi*v_0)^2*sin(i*Pi*v_0*tau)*sin(j*Pi*v_0*tau)+2*mu*(i*Pi*a_0)*sin(i*Pi*v_0*tau)*cos(j*Pi*v_0*tau)
      ):
f:=tau -> Vector
      ( n,
        (i,j)-> x*Pi^2*(Pi/2+v_0^2*sum(((1-(-1)^k)/k^3)*sin(k*Pi*v_0*tau),k=1..infinity)-a_0*sum(((1-(-1)^k)/k^4)*cos(k*Pi*v_0*tau),k=1..infinity))*sin(i*Pi*v_0*tau)
      ):
      
 
X:= Vector(n, i-> x[i](t)): 
sys:= M(tau).diff~(X, tau$2)+C(tau).diff~(X, tau)+K(tau).X=f(tau):