The pde is fourth order in X, fourth order in theta, and second order in t, which means that it needs 4+4+2=10 boundary/initial conditions for an explicit solution.

Your suggestion

sol := pdsolve(PDE1,'build');

which should be 'PDEs' not 'PDE1' returns an expression with 10 arbitrary constants.

Given tha the OP only has 6 boundary/initial conditions - I think we are on a loser here

Nice,thank- but shouldn't this be the default?

I got the "analytical" solution using the command in my original post. However this is almost certainly wrong. ThomasRichard has pointed out that you are using

x, X,  theta, t as variables: (x occurs in bcs3)

I suggest that first of all you decide which are the independent variables for the problem.

Bear in mind my first comment - that Maple can only (numerically) solve PDEs in two independent variables

You are in luck! Your pde can be solved analytically!

Try

sol1:=pdsolve(`union`(PDEs, bcs1, bcs2, bcs3));

whihc produces

sol1 := w(X, theta, t) = ((10*cos(1)+10*sin(1)-10*cos(X-1))*exp(2*x+1+X)+(10*cos(x-1)-10*cos(1)-10*sin(1))*exp(x+1+2*X)+(10*cos(X)-10)*exp(2*x+X)+(-10*cos(x)+10)*exp(2*X+x)+(10*exp(X+x)-10*exp(2+x+X))*sin(X-x)+10*sin(1)*exp(1+2*X)-10*sin(1)*exp(2*x+1)+(-10*cos(X)+10*cos(x)+10*sin(X)-10*sin(x))*exp(2+x+X)+(-20*cos(x-1)+20*cos(X-1))*exp(x+1+X)+(-10*cos(X)+10*cos(x)-10*sin(X)+10*sin(x))*exp(X+x)+(10*cos(1)-10*sin(1)-10*cos(X-1))*exp(X+1)+10*exp(x+1)*cos(x-1)+(-10*cos(1)+10*sin(1))*exp(x+1)+(10*cos(X)-10)*exp(2+X)-10*exp(2+x)*(cos(x)-1))/((-cos(1)-sin(1))*exp(2*x+1+X)+exp(2*x+X)+(-cos(x)+sin(x))*exp(2+x+X)+2*cos(x-1)*exp(x+1+X)+(-cos(x)-sin(x))*exp(X+x)+(-cos(1)+sin(1))*exp(X+1)+exp(2+X))

If you are happy with the idea that connecting the dots sometimes draws lines within the sphere, then the attached works for me

geodesic.mw

If you actually want to draw the connecting lines on the surface of the sphere, then you have to construct a sequence of spacecurves which are greatCircles (aka geodesics) between each successive pair of points. This is trickier than it looks.

I'd develop this on a new worksheet just by working on two points and trying to produce a procedure which generates the great circle between them.

In order to convert to a first-order system, Maple has to be able to isolate the highest order derivatives in each equation. The term

(diff(f(eta), `$`(eta, 3)))^n

with n=2 causes the problem, because isloating (diff(f(eta), `$`(eta, 3))) is a quadratic problem with two solutions.

One can generate both of these solutions and producce two separate systems of ODEs (only you know which one you want!)

When I try to solve either/both of the resulting ODE systems, I get a singularity error at the right end point ie eta=inf=5.

When I try to get around this by using a midpoint method, I get a datatype error because complex values are being generated for both systems.

See the attached

uglyODEs.mw

Before I go any further I'd like you to check your original equations and BCs very, very, carefully.

What did I do -well

1. Some general syntax tidying
2. Commented out the Pade approximation stuff: I think this was just making things more complicated
3. Can't produce solutions when the limit at y=infinity is used, because this limit is undefined. So start with the limit at 10¹, compute the solutions, then redo for p=2,3,...and until the solutions are "no longer changing" - ie the current list of values is within abserr (which I defined as 1e-06) of the previous set.

Solutions obtained were

[[A = -.7826356580, B = .1991533924],
 [A = -.9612853602, B = .9613743694],
 [A = -4.878194845, B = 6.777427042],
 [A = 1.390879992, B = 11.22909487],
 [A = -4.812628362+5.366978521*I, B = -2.705290847+13.44775887*I],
 [A = -.3687462685-.7543637171*I, B = -1.864599186+2.915730929*I],
 [A = .7498822399+0.8909367120e-1*I, B = -2.842643414+1.790102688*I],
 [A = 1.204969827-5.296330878*I, B = -11.38132992+.1118543744*I],
 [A = .8564784932, B = -.3590096093],
 [A = 0.6363850082e-1, B = -2.949292641],
 [A = -6.265195047, B = -70.43835456],
 [A = -2.726994841, B = -75.44978462],
 [A = -10.16962853, B = -103.3206546],
 [A = -2.505843017, B = -109.4092818],
 [A = -2.536884197, B = -435.2896991],
 [A = 1.204969827+5.296330878*I, B = -11.38132992-.1118543744*I],
 [A = .7498822399-0.8909367120e-1*I, B = -2.842643414-1.790102688*I],
 [A = -.3687462685+.7543637171*I, B = -1.864599186-2.915730929*I],
 [A = -4.812628362-5.366978521*I, B = -2.705290847-13.44775887*I]
]

See the attached

DTMprob.mw

http://www.maplesoft.com/support/help/Maple/view.aspx?path=alloc

Particularly the last paragraph.

Since all of the 'consumer' packages (ie maple/macsyma/matlab/mathematica) seem to use the same underlying routines from NAG/LAPACK/BLAS, for LinearAlgebra/Optimisation problems I'd be sceptical about any of these which claims to be "better" than any other in this field.

My scepticism would extend to GAMS/Lingo

I have read the sentence

Export puts this proc into an area that doesn't seem to show up when the document is compiled.

and I have no idea what this means!!

If I just do a quick check of the Maple help for "mapleinline", then on the help page ?Printing and Viewing Maple LaTeX Documents, It states

The LaTeX macro \mapleinline is used to encode Maple 2-D inline math displays.  For example,
  \mapleinline{inert}{2d}{sin(x^2);}{$\sin(x^2)$}
corresponds to an inert Maple expression that appears in the Maple worksheet in 2-D mode. The Maple command sin(x^2) is displayed in LaTeX by "$\sin(x^2)$".

Seems clear enough - what exactly is your problem?

The OP's original request was how to construct fractional differential equations, for which I pointed him/her to fracdiff. My only concern when doing this was: how would OP solve the resulting fractional DE system??? Being able to write the equations with no hope of solution, is rather pointless.

So I dug into Maple a little further, and discovered that it appeared to be able to generate the Laplace transform of fractional derivatives (I was surprised!). This suggested one possible solution method - perform the Laplace transform, solve the set of algebraic equations, perform the inverse Laplace transform. This is three-step process and any one of these steps could fail. However it was a possibility, so I though I'd mention it

If you have three dependent variables [h1(), h2(), h3()] then you need three differential equations

Looking at the set of ODEs in dsolve2.mw, it is 4th order in h1(), 3rd order in h2(), and 4th order in h3() - so you need 4+3+4=11 boundary conditions

In dsolve2.mw, you have three ODEs and eight boundary conditions: this cannot be solved, because you have an insufficient number of bcs

In your original worksheet dsolve.mw, you end up with two ODEs and 11 boundary conditions: this cannot be solved because you have an insufficient number of ODEs (which I should have spotted in my original response - sorry!)

Well  it is possible to create a fractional differential equation using fracdiff(), but could one solve it!???

I'm a little out of my comfort zone here, but I know that the simplest(?) solution method is to use the "fractional" Laplace transform. Somewhat to my surprise, Maple seems to produce the correct answer, for these, since

inttrans:-laplace(fracdiff(y(x),x,1/2),x,s);

actually works!

So I guess this mught be a viable approach. Good Luck!

In order to numerically solve a system of higher order ODE's, the first stage is to convert the system to one involving only first order ODEs, accompanied by the introduction of a lot of new variables - see any decent textbook on ODEs or even the wikipedia article https://en.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations which clearly states

Without loss of generality to higher-order systems, we restrict ourselves to first-order differential equations, because a higher-order ODE can be converted into a larger system of first-order equations by introducing extra variables. For example, the second-order equation y'' = −y can be rewritten as two first-order equations: y' = z and z' = −y.

Maple allows you to examine this stage in the solution process using the DETools:-convertsys() command, which ought to return a (augmented) system of first-order equations. However

DETools:-convertsys(SYS2, bcs2 union {bcs22}, [h1(theta),h2(theta),h3(theta)], theta);

eventually (after about ~10mins on my machine) returns the error message

Error, (in DEtools/convertsys) invalid specification of initial conditions

Most of your boundary conditions are pretty standard - so I suspect the problem is with bcs22, which evaluates to

-.35329*(D@@2)(h2)(1)-0.16526e-4*(D@@3)(h1)(1)-0.17833e-2*(D@@3)(h3)(1) = 0

which isn't so much a boundary condition as a "combination" of boundary conditions It is not obvious (to me) why this bc should cause an issue, but I strongly suspect that it is the source of your problem.

Any chance that you could come with an alternative bc?

As a starting point I decided to evaluate yur requirement with the assumption that all multiplications are commutative using

restart;
eq1:=x*y-q*y*x=1;
f(x,y)=simplify(subs(isolate(eq1,x), x*y^n-q^n*y^n*x));

Somewhat to my surprise, this produced the answer you desire. I find it unlikely (but I guess possible) that commutative and non-commutative multiplications would result in the same answer - so are you really sure about the non-commutative requirement??

Maple does have a non-commutative multiplication operator, '.', but this is primarily intended for matrix/vector products and, as the help page says

Arguments that are not of type Matrix, Vector, constant, or Array are ignored, and A.B remains unevaluated.

Hence I suspect that Maple would be unable to produce the output you want, with the restriction of non-commutative multiplication :-(