@Subhan331 I ask again: Why do you need Maple code for an algorithm that produces nonsense?  I wouldn't waste my time on that, and nor should you.

As a student and a beginner, you should spend your time on more productive pursuits.

I have no problem with "with(combinat)" in Maple2017 and Maple2020 on Linux.  I suspect that you have typed something incorrectly.  Post your worksheet so that people can examine the source of the error.

To post a worksheet, reply to this message, and note the big fat green arrow in the toolbar of the window where you compose your reply.  Click on that arrow to upload your worksheet.

@Subhan331 Doing the exercise that I noted takes just a few lines of simple calculations by hand.  What's the point of writing Maple code if the algorithm produces nonsense?

As to fractional order derivatives, isn't is too soon to jump into that?  See if you can make sense of the algorithm applied to traditional PDEs first.

An alternative to dharr's suggestion is to press Ctrl+Shift+G (that is, the Ctrl, Shift, and G keys together) followed by the Greek letter code.  For instance, Ctrl+Shift+G followed by e produces ε, and Ctrl+Shift+G followed by s produces σ.  Most Greek letters map to their corresponding Latin keys.  For the complete table of Greek character, look up the help page 2-D Math Shortcut Keys and Hints.

The Ctrl+Shift+G combination works on Linux and Windows.  On a Mac, type Command+Shift+G instead.

@Preben Alsholm What you have quoted can be an explanation of what we are seeing.  It's not quite clear to me how, but it seems relevant.  Thanks for your input.  Because of that I have converted your Reply to an Answer.

@Pascal4QM You need to say
A,C := test()[[1,3]];  

@mmcdara  Oh, I see what you mean—you want to consider the case where the beam is not mechanically attached to a support, and therefore it may detach and lift off.

That requires handling boundary condition given as inequalities rather than equalities.  That's a difficult problem, both mathematically and computationally.

However, there is some hope.  Beamsolve reduces Euler's PDE to a system of ODEs.  It is conceivable that the one-sided constraints at the supports may be handled through dsolve's Events mechanism as described in the help pages:

?dsolve,events
?How Do I Solve an Ordinary Differential Equation?

That can be a fun project for someone to take up.
 

@mmcdara Any of the values of u, u_x, u_xx, and u_xxx may be specified for the boundary conditions.  I have picked u=0 and u_xx=0 as defaults because that corresponds to a simple support (a ball-and-socket support, as you put it) which I assume would be the most common type of support. I don't understand what you mean by "this prevents the simulation of simple supports", since the user's specification overrides the defaults.  Let me know if I have misunderstood you.

@rookie Beamsolve computes the values of u(x,t) in order to plot the solution but it does not provide a user-accessible function to extract those values.  I will think about adding one but this may take some time.  In fact, it will be good to have functions that return not only u, but also u_x, u_xx, and u_xxx.

@mmcdara Thanks.  Calculating the stress is easy -- it's just a multiple of u_xx.  Producing good-looking graphics to highlight it, is another story.  That can wait until after I am retired :-)

@Pemudahijrah01 Okay, then try this:
DEtools:-DEplot3d(sys, [x(t),y1(t),y2(t)], t=0..500, stepsize=0.5,
    [[x(0)=0.6, y1(0)=0.3, y2(0)=0.3]], arrows=cheap, linecolor=black);

@Pemudahijrah01 Is this good?
DEtools:-DEplot(sys, [x(t),y1(t),y2(t)], t=0..50, stepsize=0.1,
    [[x(0)=0.6, y1(0)=0.3, y2(0)=0.3]], scene=[t,x(t)], x=0..0.6,
    linecolor=black, thickness=1);

@Pemudahijrah01 When you change Mu[2]=0.74 to Mu[2]=0.38, the equilibrium shifts to  {x = 0.5233951682, y1 = 0.2697791086, y2 = 0.}.  Note the y2=0 part!

The new equilibrium is still stable.  The x(t) and y1(t) converge quickly to their equilibrium values, but y2(t) takes a very long time to go there.

To adjust the phase portrait for the new parameter, change the previous
y2min, y2max := 0.1, 0.2;
to
y2min, y2max := 0.0, 0.01;
Then we get:

It helps to plot x(t), y1(t), y2(t) as functions of time.  These show that x(t) and y1(t) stabilize within 0 < t < 50, but for y2(t) you need 0 < t < 1000.

DEtools:-DEplot(sys, [x(t),y1(t),y2(t)], t=0..50, stepsize=0.1,
    [[x(0)=xmax, y1(0)=y1max, y2(0)=y2max/2]],
    x=xmin..xmax, y1=y1min..y1max, y2=y2min..y2max, scene=[t,x(t)],
    linecolor=black, thickness=1);

@rookie Here is a version of the worksheet in which I have replaced some of the newer Maple syntax with their older equivalents.  This one works with Maple 2017 (tested) and perhaps earlier (not tested).

Download worksheet: euler-beam-with-method-of-lines-2017.mw