It's not obvious which variables to eliminate, but this works. 

Download elim.mw

For a numerical solution the limit boundary condition can't be used. If you want to approximate infinity by by a large number you can use, say, U[2,n](20)=0. But then you have boundary conditions at -Pi, 0 and 20; the solver needs just two boundary locations, so I replaced it with a boundary condition at 0, which you will need to modify to what you want. You need also to replace x[01] with X[01] to avoid confusion with the simple variable x. Then it is possible to get a solution.

Download Thesis_(1).mw

So here is a workaround,  which just loops as @sursumCorda suggested.

Download PrincipalPart.mw

Not sure how you want to label them, but here are a couple of simple possibilities, using the caption as a label.

Download Plotlabels.mw

@lcz For IsSubgraphIsomorphic, Maple uses a constraint algoritham and SAT solver, perhaps this algorithm? doi 10.1007/s10601-009-9074-3 or here? doi: 10.1016/j.artint.2010.05.002, but I'm not sure there is an equivalent for the induced case.

The VF2 algorithm and its (much) improved variants VF2plus, VF2++ and VF3, seem to be widely used. Here I only tried VF2, and didn't try to optimize the data structures (mainly sets straight from the paper); probably moving to one of the improved algorithms would be the next step. I didn't check it on a large number of cases; perhaps you have some other cases to test it.

[Edit2: Updated version here now implements some parts of the VF2++ algorithm, and removes redundancy in search tree]

It seems fast to find matches if there are some, but of course is slower to show there are no matches. Hope this is useful.

[Edit - full VF2++ below]

Download VF2conndegAlgorithm4.mw

I changed your code to

f := "this:///Images/Maple.jpg";
img := Read(f);

and it works - you then get some warnings later that rotation angles should be in radians, which I'm sure you can fix.

Although

diff(ln(GAMMA(x)), x)=Psi(x)

Using the chain rule we find

diff(ln(GAMMA(1/x)), x)=-Psi(1/x)/x^2;

which explains the missing -x^2.

I don't think there is a builtin command, but there an implementation of a heap, which allows it to be done easily. If you wanted the whole list partially sorted, then just sort the selected ones and follow with the rest. (The smallest ones are in decreasing order; the largest are in increasing order.)

Download partselect.mw

One solution to this error is to supply an approximate solution, and since you said tanh(x) was a known solution, I tried that. But then I realized tanh(x) goes to -1, not 0, as x->-infinity. If I change the boundary condition to z(-15)=-1 it works. But if you really wanted z(-infinity)=0, then you can try a better approximate solution.

dsolve.mw

See

https://www.mapleprimes.com/questions/235168-How-Do-I-Generate-Magic-And-Semi-Magic

for some solutions.

You can add a constant to each cell to get other ones, but not sure what exactly you mean by random.

For your first case, the initial conditions are specified as 

dsolve({DE, R(0) = 1, D(R)(0) = 1}, numeric, range = 0 .. 20)

I made up a value for the derivative in the second condition; you will no doubt have a better value. Or perhaps you wanted a boundary condition as your second condition.

And a similar problem applies for your second problem.

PS: In your second problem in defining F you probably wanted an explicit multiplication after the first ), so (...)*(...)

Your second problem doesn't seem well-posed at theta=0.

Here some progress:

DIFFERENTIAL_EQUATION.mw

Download pade2.mw

The degree 3 result here is also a problem - I'll submit an SCR.

Not as general as @acer's solutions. It seems evalhf handles the Gamma function but not factorial so in this particular case converting to GAMMA works.

[Edit: This explanation is incorrect - see below]

failed_plot.mw

[Edit: the OP has now completely changed this post, so my responses here and below are no longer relevant.]

The equations in your post outline a laplace transform analytical solution and you seem to want to follow that, but in your worksheet you set up for a numerical solution of the original untransformed equations, abandoning the laplace transform method.

So forgetting laplace transforms and following your worksheet, you have two errors in formulating a consistent problem. You had strange parentheses after alpha0 in your Cond, which I guessed as a correction (in red). [Edit] What about initial conditions? 

new_paper_dust.mw

If you want the analytical solution then you can start by transforming the pdes to odes with intrans:-laplace, gving you Eq (16) etc, then solve the ode system - but then it is not a good idea to give numerical values to your parameters.

Maple tries to figure out the file format you want from the file extension, but doesn't recognize .dat. You can use the format option to specify which format you want, for example CSV:

Download testfile.mw

Edit: for saving Maple objects in a .m file for reading in later into Maple use the save command