You gave h in terms of a in the line before the plot, so now f is an expression in x and a, and you need a range for a. So change h to a in implicitplot and you get a plot. 

In general to diagnose these sorts of errors, just change the command to an invalid one (I changed implicitplot to iimplicitplot) and you can see what the command actually received.

Well, my answer is more in the context of fitting data to equations, but it could be used here. Notice that the distance between two points z1=x1+I*y1 and z2=x2+I*y2  in the conmplex plane is abs(z1-z2) = sqrt((x1-x2)^2+(y1-y2)^2) so the square of this is abs(z1-z2)^2 = (x1-x2)^2+(y1-y2)^2. So if you add these for all your data point pairs you get an overall sum of squares of the real and imaginary deviations, which is a measure of the overall agreement. (This is what you would minimize in a comnplex least squares fit)

As @mmcdara has shown, MathML is how to do it. But to find the MathML , you can create what you want in 2D math and then convert to 1D math. For example, in 2D math use w^2 to create a superscripted w. Then select the w^2 and use the right-click menu (ctrl-click on Mac) to do 2-D math -> convert to -> atomic variable. Now use 2-D math -> convert to -> 1-D math input, and you hve your variable:

`#msup(mi("w"),mn("2"))`

For combined sub and superscript, use the Layout palette to generate it in 2D,

The SparseIterative method requres shape = symmetric, storage = sparse, datatype = float[8] and is for real symmetric or complex Hermitian Matrix (so the Matrix you enetred was changed to be symmetric)

Sparse.mw

Edit: methodoptions = [itermethod = CG, precon = false]  leads to the symmetric requirement (error message without) - in fact CG is for positive definite.

Your immediate problem is the Q's in the denominator. If you take the numerator and set it to zero, then you can almost get the answers you want PROVIDED P=1, which didn't seem to be in the problem statement. And why is there a d_1?

coeffs.mw

Maple has not found a solution, but has partly solved the problem (in the second solution with the _a and _b(_a)) by reducing the system of two ODEs down to one. Where you go from here depends on whether you are OK with a particular solution or not, or if you expect a solution in a particular form (e.g. a polynomial). You also have an unknown _C1 which you can eliminate by providing an initial/boundary condition (which might also simplifiy things).. To build a solution from the partial one, you can use buildsol in the DEtools package - see its help page.

As you probably know, you should be able to get a numerical solution even if you can't get an analytical one using dsolve(...,numeric)

LibraryTools:-ConvertVersion converts between .lib and .mla - see the help page.
 

Maple doesn't seem to know it is separable in my version (2017), but the steps you suggested can be carried out and Maple can solve the ODE.

PDE.mw

Allthough there are no explicit solutions, you can still use allvalues to plot them (only first two are real) or process them further

RootOfQuestion.mw.

The procedure exits through the return statement (giving the 0.5 output), so code after that is not executed.

ps. array should be Array for moderm Maple.

@DarkMath Your examples

f(h(r)) ==> int(h(r),r)
f(h(r+a)) ==> int(h(r+a),r)

require different code. The second one somehow has to know that "r" and "a" have different roles, otherwise you will get

int(h(r+a),r+a)

Your first one can be done as

f:=proc(g::function)
   int(g,op(1,g));
   end proc;

then f(h(y)) gives int(h(y),y)

But then you can't pass an expression to f, so not sure if this is exactly what you want.

Edit: For your second one, if it will always be "a", you can do

f:=proc(g::function)
   int(eval(g,op(1,g)=op(1,g)+a),op(1,g));
end proc;

Then f(h(y)) gives int(h(y+a),y)

I just downloaded the latest MaplePlayer (Windows 64 bit, 2020.1) and loaded a Maple Worksheet with a text area component, and was able to select the contents and use ctrl-C to copy. It pasted into a text file (empty notepad file) as I expected. So perhaps this is a version issue or you wanted to copy other types of content?

use 

infolevel[sum]:=5;

Then enter the sum (rather than Sum) and some information is given:

sum: process the input arguments
sum: definite sum
ClosedForm: try the front-end to compute simple definite sums
ClosedForm: simple definite summations
ClosedForm: try the method of integral representation
ClosedForm: try the method of first computing the indefinite sum
ClosedForm: infinite summation
Indefinite: process the input arguments
ClosedForm: try the method of converting the definite sum to hypergeometric functions
ClosedForm: definite summation using hypergeometric functions
ClosedForm: definite summation using hypergeometric functions

So Maple has some routine that does definite summation using hypergeometric functions, which is what you deduced already I believe. You can trace the actual calls and figure out more using "trace", but following the Maple code can be difficult.
 

dsolve directly gives the general solution in terms of unspecified constants _C1 etc. For each initial condition you specify there will be one fewer unknown constant. It's not clear to me exactly what else you want.

GeneralSolution.mw

Well, I don't think it should hang. On the other hand if you try it (Maple 2017) without assumptions then it gives GAMMA(n+1) (I didn't need method =MeijerG). And with an extra step you can get Maple to tell you what you need, a bit like Mathematica.

FunctionAdvisor(GAMMA); under definition shows that GAMMA(z) is defined for Re(z)>0, and under singularities shows:

Or(z::nonposint, z = infinity+infinity*I)

So this can be a guide to what assumptions would work.