Your program refers to PlyTable but as far as I can see, it cannot read its value.  To verify that, print the value of the CurrentAngle within your main loop, as in

for i from 1 to NumberPlies do:
    CurrentAngle := PlyTable1[i,1]:
    printf("i = %d, CurrentAngle = %a\n", i, CurrentAngle);
    ...

 By the way, why do you have a colon after "do"?  Remove it.

A general advice: There is an old saying to the effect that for any problem that you cannot solve, there is a simpler problem that you cannot solve.  Try that one first.

Your current program is too complex to debug. I suspect that you received it from someone else and are trying to modify it.  If that's the case, you will be better off by writing your own version, from the beginning and without bells and whistles.  Start simple.  Do away with widgets that can hide data.  You may make your program fancier later, after you have established a foothold and verified that the plain version works.

@acer The solution is to remove the f(x) := whatever syntax from Maple altogether.  As a help to new users, an attempt to f(x) := whatever should trigger an error message with a text like "did you mean f := x -> whatever?".

@acer Thanks for the explanation.  That's indeed an unexpected behavior (a polite way of saying a "bug").  I still don't understand why gamma(1) evaluates to −0.07281584548.  Do you know where that number comes from?

@dearcia Maple is unable to plot the arrows in your original version because the vector field is undefined on the planes x=0, y=0, and u=0.  To avoid the singularity, don't go near those planes. Adjust the plotting region like this:

 x = 0.1 .. 1, y = 0.1 .. 1, u = -1 .. -0.1

@PhD_Wallyson I don't have access to articles 1 and 3.  Send them to me and I will have a look.

@ Oh, I see, this model is a lot more complicated than I would have expected.  I wouldn't be able to devote the necessary time to go into its details.  I hope that someone else in this forum will.

@Paulo Baumbach  That's a good description of the problem and I have a better understanding of what you are attempting to do, but some details are still fuzzy.  For instance, I can't tell whether "phase 1" and "phase 2" refer to intervals in time or in space.  Also I don't see the differential equations.  I don't expect you to explain all that here because that would amount to copying  a good deal of the article you have referred to.

It is a hoiiday weekend in the United State and the university library is closed.  I will see if I can find the article next week and see what it is about.

I looked through your lengthy worksheet but was unable to understand the main idea, and therefore I am unable to offer a suggestion.

It will help if you could provide a statement of the problem in words and mathematical equations, separating them from the computational details.

@ahmadtalaei What is the purpose of introducing the variable eta into the operator Do?  I cannot follow that line of thought.  See if you can clarify what it is that you want to do.

As to D, D^2, D^3,  they are the equivalent of the notation w', w'', w''' in calculus.  Thus, the ODE that I have obtained looks like this:

A(η)*w''''(η) + B(η)*w'''(η) + C(η)*w'(η) = 0

After you have calculated the solution w(η), you may plug in r*sin(θ) for η. to obtain w(r*sin(θ)).  There is no differentiation involved in that operation, and therefore the question of "when taking the derivatives of " should not arise.

@Thomas Richard I didn't know about that MapleSim animation. Looks good!  Thanks for pointing it out.

@Preben Alsholm Oops, I was careless.  I have corrected the name now and thus have baptized him Dutch.  Sorry ;-)

@PhD_Wallyson Have a look at:

https://www.mapleprimes.com/posts/213083-Eigenfunctions-Of-A-Multispan-EulerBernoulli?sp=213083

See if that helps you with your work.

@PhD_Wallyson I quickly looked over your worksheet.  Perhaps it's possible to fix it to make it work but altogether the approach there is somewhat disorganized and not very pleasant to work with.

I will do a writeup on how to apply the Krylov-Duncan functions in a systematic way to analyze multi-span beams, and will let you know when it is ready.  I may take a day or two.

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