The idea of "the arithmetic mean" in this setting is ill-defined.  Let t be the angle that the radius makes with respect to the x axis, and let s be the arc length measured along the ellipse.

The length L of the radius may be expressed as a function of t or s.  If L is viewed as a function of t, for example, then the arithmetic mean would be 1/(t2-t1)*int(L(t), t=t1..t2).  On the other hand, if L is viewed as a function or s, then the arithmetic mean would be 1/(s2-s1)*int(L(s), s=s1..s2), which is something quite different.

L also may be viewed as a function of x, or y, where (x,y) is the coordinates of the tip of the radius.  You can think of other variables that can serve to parametrize L.  There is a different "arithmetic mean" corresponding to each choice.

Although you don't state this explicitly in your worksheet, it is clear that the intention is to solve the partial differential equation of a vibrating Euler's beam through separation of variables.  It is also clear that ω² is the separation constant that arises in the process.

Things go wrong when you assign a value to ω², that is, where you set ω²=a² c. Where does that come from?

To fix, go over the process of separating the variables carefully.  Additionally, examine the physical assumptions that lead to the mathematical model.  Specifically, the motion of the beam certainly depends on how you get the motion started.  At no place in your worksheet you have said how the motion gets started.  Look up the term "initial condition" in your reference book.

@jefryyhalim I don't know how to compute the next eigenvalue but I doubt that it would be useful.  Generalizing from the 1-D model, the higher eigenvalues and their associated buckled shapes are likely to correspond to unstable equilibria and therefore you would not see them in experiments.

If your actual experiments show wrinkled shapes, then perhaps your mathematical model is at fault.  As I have said earlier, a plate is inherently a 2-D object and the proper model for its deformation is a partial differential equation (PDE).  I have no idea how you go from there to your system of ordinary differential equations (ODEs).  If I were to investigate the source of the mismatch between computations and experiment, I would start with a close examination of the derivation of the mathematical model.

@jefryyhalim As Carl has correctly pointed out, my suggestion for the midway-point condition V1p(L/2)=1 is not tenable.  (Such a condition works with symbolic solutions, but not with numeric.)  An alternative is to specify a condition at an end point, such as D(V1p)(0)=1.  You need to decide which of the dependent varaibles is suited for such an assumption.  The numerical value (=1 in this case) is immaterial as long as it is nonzero.  You will get the same eigenvalue sigma regardless of that value.

The following simple illustration may be of some help.  The numerical solver picks the a=1 eigenvalue.  (Actually this problem has infinitely many eigenvalues.)

Download zz.mw

@jefryyhalim I see what you are attempting to do.  In your worksheet Longitudinal_ODE.mw you have a system of 5 fourth order differential equations and 20 boundary conditions.  All of these are linear and homogeneous.  Clearly there is a trivial solution consisting of all zeros.  (This was pointed out by mmcdara.)  Furthermore, from your original worksheet SP1JH1B1.mw I gather that you wish to introduce a parameter, sigma, into the coefficients of your system, and then you wish to determine sigma so that the system has a nontrivial solution.

I don't know the physical meanings of your variables Vlp1, Vup1, etc., but you know what they mean, and probably you can guess what a buckled configuration should look like.  Let's say that you expect that the solution Vp1(x) [you pick the right one] to be something like an arch, that is, zero at the endpoints and nonzero in between.  Then add an extra condition, V1p(L/2)=1 to your existing 20 conditions.  Then solve your system of differential equations together with the 21 conditions numerically, for the six unknowns consisting of your original five, plus the sigma, as in:
dsol := dsolve(
        { the five DEs } union {the 21 conditions},
        { Vup1(x),Vp1(x),Vp2(x),Vlp1(x),Vlp2(x), sigma },
        numeric, output=listprocedure);

If this succeeds, then sigma will be the desired buckling load.

Good luck!

@jefryyhalim I don't quite understand your mathematical model.  The buckling and deformation of a plate are customarily modeled as partial differential equations with x and y coordinates denoting the position on the plate.  Your equations are ordinary differential equations in terms of a single independent variable x.  Where is the plate's y coordinate in your model?

Question: To analyze the buckling, you need to specify boundary conditions all around the plate.  That involves both the x and y coordinates. What are the plate's boundary conditions?

Suggestion 1:  Before analyzing the buckling of  a plate, I suggest that you do the much simpler problem of buckling of a beam. Do you know how to do that?

Suggestion 2: The very last command in your Longitudinal_ODE.mw is:

sol_L:=evalf(value(dsolve(...)));

That composition of three commands obscures the ensuing error message since it is unclear which of the three commands has led to it.  Don't compound commands. Execute them one at a time to ensure the success of one before applying the next.  Thus, we begin with the innermost command:

dsolve(...);

and we immediately see that it returns nothing, indicating that there is something wrong with your system of differential equations.

@jefryyhalim Here is how one solves a differential equation along with initial conditions in Maple:

de := diff(y(x),x,x) + 9*y(x) = 0;
ic := y(0)=1, D(y)(0)=2;
dsolve({de,ic});

This yields:

This is illustrated in the very first example of Maple's help page on dsolve(). You should read that page.

You say you need the coefficients in order to form a matrix and then calculate eigenvalues.  That strikes me as a roundabout way of achieving something that may be doable in a shorter way.  If you explain your overall goal (rather than the steps toward it) someone may be able to offer constructive advice.

Another remark: In several places in your worksheet you have constructs similar to

f(x) := x^2;

That doesn't do what you want.  To see the problem, do f(7).  Do you get 49?  No!

To correct way is:

f := x -> x^2;

Try f(7) now.

I looked over your worksheet but couldn't make heads or tails of it.  If possible, just write down the complete set of your differential equations and boundary conditions clearly as plain mathematics (no Maple calculation).  Then someone here should be able to help you code that in Maple.

One observation: In your approach to solving a differential equation, you obtain the general solution first, and then apply the initial/boundary conditions in order to determine the integration constants.  Maple is capable of doing those two steps together.  Doing so will simplify your calculations significantly.

@Carl Love What you have written is correct, however it doesn't address the original poster's question who is asking why the limit of the difference quotient yields zero rather than 3.  I don't have a good explanation.  Is it perhaps due to a bug in Maple's limit()?

What you have shown looks good.  I suspect that the error is due to something else that you have not shown.

Consider posting your entire worksheet.  For that, click on the big fat green arrow which appears in the window where you type your message to MaplePrimes.

@verdin OK, I will do the calculations and post it as an Answer to close this thread.

@verdin Have a look at this modified worksheet where I have picked arbitrary parameter values and attempted to find the system's equilibria.  No equilibria are found.  That's because the parameter values have been picked blindly.  With your insight into the origins of the problem you should be able to come up with a meaningful set of parameters which might yield an equilibrium.

mw.mw

@tomleslie The OP did not ask for a path integral.  He wants to integrate the given expression over (p,q,w) ∈ R³×R³×R³, where p=⟨p₁,p₂,p₂⟩, etc., are independent variables. Thus, we have a nine-fold integral over R⁹.  I haven't attempted doing that but just from the looks of it I expect it to be complicated.

@tdavid I have fixed a few errors in your worksheet.  Here is what we get now.

Download mw.mw

@verdin It appears that the steady-states of your system are expressed through the roots of a complicated transcendental equation.  I am afraid that there are no analytic expressions for the roots.  You may have to analyze stability for specific numerical choices of the problem's parameters.