@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

@Neel See if this helps.

Download eigenfunction-calculation.mw