You have a system of 15 first order PDEs in 15 unknowns, each being a function of z and t.  You discretize the z variable and obtain a system of 15n ODEs in 15n unknowns.  Even for moderate n that's a very inefficient way of solving the system, as well as being a very poor approximation.

 I suggest applying pdsolve() to solve the system of PDEs directly.  Have you tried that?

@Neel There are many different approaches to calculating the eigenfunctions.  I showed the simplest approach in my earlier reply.  Since you are asking for a solution with with matrices, here is the matrix version.  As to the forced vibration, I showed how to to that in my pdf writeup.

Download eigenfunction-calculation2.mw

@Neel See if this helps.

Download eigenfunction-calculation.mw

@Neel See if the writeup in the attached PDF answers your questions.

By the way, in your line of work you will need to go beyond writing mathematics by hand.  Look into learning LaTeX if you haven't done so already.  My writeup is done in LaTeX.

steady-state-oscillations.pdf

Note added later:
I spotted a few (harmless) typos in the previous PDF.  Here is one with corrections.

steady-state-oscillations2.pdf

You are suffering from culture shock.  You must already know that when in a foreign country, speaking louder in your own language does not make you any better understood.  Learn the local language and customs instead of complaining that they do not agree with yours.

@Neel In equation (9.B.20) you have the force expressed as a delta function.  The picture below that equation also indicates that you do mean a delta function.  But now you are saying q(x,t)= Fo*(1-(x/L))* cos(omega*t) without a delta function.  You need to decide which one you want.

@Neel I looked quickly through your PDF.  The calculations up to (9.B.19) appear to be correct, but after that point I felt lost.

In (9.B.22) we have the equation of motion:

ρ A w_tt + E I w_xxxx = q(x,t)

which is correct, but we have no clear statement of what q(x,t) is.  You need to say
q(x,t) = something.

What is that something?

The boundary conditions
w(0,t) = w_xx(0,t) = w(L,t) = w_xx(L,t) = 0.
for pinned ends are good.

I also don't see what the initial condition are.  You need to specify the initial displacement and initial velocity as

w(x,0) = something
w_t(x,0) = something

Can you tell me what those "something"s are?

@Earl   The Harvard Weekly Problems is a good collection math and physics puzzles.  The one you are referring to is Problem #15.  Beware though that the shape discussed in that problem is the 3D solid obtained by rotating your lamina about the horizontal axis.

@acer  As you have observed, computing  int(int(x,y=-f(x)..f(x)),x=0..1) without the method=meijerg option produces an answer in terms of EllipticE and EllipticK functions with imaginary arguments:

int(int(x, y=-f(x)..f(x)), x=0..1);      
                                     
                      -4/5 EllipticK(I) + 4/5 EllipticE(I)

As we have seen in my original replay, that answer then gets reduced to a real value after applying convert and simplify.

Interestingly, if we change to polar coordinates by letting x=r*cos(t), y=r*sin(t), the equation of the lamina's boundary changes to r = sqrt(cos(t)), where t goes from −π/2 to π/2, and the integration produces an answer in terms of EllipticE and EllipticK functions with real arguments:

int(int(r*r*cos(t), r=0..sqrt(cos(t))), t=-Pi/2..Pi/2);                       
                                   1/2                        1/2
                    1/2           2            1/2           2
              -2/5 2    EllipticK(----) + 4/5 2    EllipticE(----)
                                   2                          2

simplify(convert(%, hypergeom));                                              
                                   1/2           2
                                  2    GAMMA(3/4)
                              2/5 ----------------
                                         1/2
                                       Pi

That is not as good as the direct answer obtained with the method=meijerg option, but at least it has no imaginary parts.

@mmcdara That's an interesting approach.  I wouldn't have thought of applying statistical tools to solve a mechanics problem.

I have converted your "Reply" to an "Answer".  Thumb up!

@Joe Riel That's good. Thanks!

@vv If we could distribute ln over a product, it would convert it to a sum and then the thing will simplify.  But I don't know how to distribute ln even in the simplest case of ln(b/a).  How does one get ln(b) - ln(a) out of this?

@PhD_Wallyson I have made a few corrections to the worksheet and now you can see the modal shapes.  Your interface conditions appear to be correct but I have not examined them carefully.  Be sure to check!  Also be sure that the values of the constants are provided in consistent units.

Download: corrections.mw

@PhD_Wallyson It is possible, in principle, to program Maple to produce a table like that, but I don't think that's worth the effort.  It would be quicker to run your worksheet 10 times with the desired values of the lengths and copy the results into your thesis.

@PhD_Wallyson In this latest worksheet, you handle the symbolic coefficients correctly.  That's good.  Continue doing it that way.

The interface conditions, however, don't look correct.  You have:

(D@@2)(X[1])(L[1]) = -(kr)*(D(X[1])(L[1]))+(D@@2)(X[2])(L[1]),     # equilibrium condition bending moment at x=L[1]

A beam's bending moment is EI X''.  Shoudn't the interface condition have EI in it?

One more thing.  In the worksheet image that you posted today, we see that the μ values are calculated incorrectly.  We even have a negative μ!  I expect that once you correct the missing EI issue, that problem will go away.  And then the plotting of the mode shape will work as before.

Yet another comment: The legnths L1, L2, L3 in the worksheet are in centimeters.  Be sure that kt, kr, E, I, A, etc., are in compatible units.