@AHSAN Your bc says q(y) = 0.  Your bc1 says q(y) = j .  Do you want both?

@Fancypants Yes, equation (20) looks good.

@Fancypants Okay, that's better.  Now you have an ODE for S[y](t).  The ODE involves phi(t), psi(t), theta(t), f(t), etc., with various subscripts.  Are those known/given functions?  If so, then your equation (13) is a simple second order linear ODE which is trivial to solve.  On the other hand, if you expect to determine the functions phi(t), psi(t), theta(t), f(t), etc,, then you need to provide additional equations.  In general, you need the number of equations to be equal to the numbers of unknowns.

@Fancypants In the worksheet you say that you expect partial_py1 to have terms involving the derivatives of S[y](t) with respect to py1.  But Maple hasn't been told that S[y](t) depends on py1, therefore it sets all those derivatives to zero.

There may be other errors in the worksheet but I don't quite see the aim of your calculations, and therefore I cannot suggest a remedy.

@spinoza Are you going to pass the solutions that you have obtained here as your own?  Or are you going to tell your teacher that Carl/Kitonum solved them for you?

@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.