@dharr If one wants a consistent aspect ratio other than the one given by scaling= constrained, it might be possible with the size option. For example, plot(..., size= [800, 500]) might preserve an aspect ratio of 800/500 (horizontal/vertical). I don't know how to test this, but presumably you do.

Even a Google search for that package turns up nothing.

@C_R That's a good point that the number of gridlines displayed (the radial spokes in this case) may be less than the number used for computation. I don't know if the reduction factor is documented, but requesting 49 will display (49-1)/2 = 24. ((I can see that you know this; I'm pointing it out for other readers): We use an odd number (49, 25, etc.) in the grid specification because of "fenceposting": That is, the number of real intervals is 1 less than the number of evaluation points, which are the endpoints of the intervals.) The extra evaluation points round out that spider-web effect. 

I've also given the featured surface a bit of transparency with transparency= 0.15 to allow the underlying polar grid to show through. 

R:= 2:  #max r of featured surface 
spokes:= 24:  #number of radial gridlines in base

plots:-display(
    plot3d(  #the featured surface
        [r, theta, (csc(theta)/r)^2], r= 0..R, theta= -Pi..Pi,
        coords= cylindrical, style= surface, shading= zhue, transparency= 0.15
    ),
    plot3d(  #the base
        [r, theta, 0], r= 0..3/2*R, theta= -Pi..Pi, coords= cylindrical,
        grid= [7, 2*spokes+1], color= COLOR(RGB, 0.97$3), thickness= 2.5, glossiness= 0
    ),
    plots:-textplot3d(  #the spoke labels
        [seq](
            [7/4*R, i, 0, i],
            i= (2*Pi/spokes)*~([$0..spokes-1] -~ iquo(spokes,2))
        ), 
        coords= cylindrical, font= ["times", 10]
    ),
    view= [map(`*`, -1..1, 2*R)$2, 0..20], axes= framed, 
    orientation= [40,40,0], labels= [r$2, z], labelfont= ["times", 16], 
    axis[1,2]= [tickmarks= [-R, R]]
);

@mmcdara Vote up. Although I think that the OP still needs to provide more detail, you've managed to make a remarkable elaboration upon a minimally informative Question.

Obviously, you need to provide more context. Do you have any Maple code, or perhaps even a journal article, to accompany this request?

What happens if you use Array instrad of array?

@Art Kalb What kind of optimization do you expect for a procedure that only does abs of a complex argument? I don't think any optimization is possible. 

@C_R You can control the polar "base" grid size like this:

plots:-display(
    plot3d(
        [r, theta, (csc(theta)/r)^2], r= 0..2, theta= -Pi..Pi,
        coords= cylindrical, style= surface, shading= zhue
    ),
    plot3d(
        [r, theta, 0], r= 0..3, theta= -Pi..Pi, coords= cylindrical,
        grid= [6, 33], color= white, thickness= 2, glossiness= 0,
        shading= none
    ),
    view= [-3..3, -3..3, 0..99], axes= none
);

@C_R What you're calling the "z" axis is showing what you call "theta", and what you're calling "r" is the angular coordinate. This is more clear if you look at the axes with numeric tickmarks. So, the "r" and "theta" in this plot are not the radial and angular coordinates, and thus they can't correspond with the r and theta as used in the Question (despite the OP's apparent approval of this plot). In other words, taking the variables to mean what they usually mean in cylindrical coordinates, you've plotted
r = 1/theta^2/sin(z)^2.

If you plot in cylindrical coordinates without giving the three coordinates parametrically, then the 1st argument to plot3d is r, the 2nd is the range of theta, and the 3rd is the range of z. This is determined by the argument order; the actual variable names that you use are irrelevant. This seems weird to me, but that's the way that it's been for as long as I can remember. Like you, I find it more natural to think of z as function of r and theta.

@acer Yes, and I just came back to correct it. (Just woke up, still lying in bed).

@Newie I don't have an idea for a solution yet. This Reply is just to let you know that given your new information about the symbolic matrix entries, I doubt that the eigenvector technique that I mentioned earlier is computationally feasible. However, it's still theoretically valid, which is possibly worth something. 

@sursumCorda Sorry for the mix-up. Try this one:

IsRectangular:= proc(L::anything, max::posint:= infinity)
local d, Op:= op@{op}, Nops:= nops, Op0:= curry(op, 0), n;
    if not L::':-list' then return false, 0 fi;
    for d to max do until 
        (n:= nops({(Nops:= Op@eval(Nops)~)}(L)) <> 1)
        or (Op0:= Op@eval(Op0)~)(L) <> ':-list'       
    ;
    d >= max or max=infinity and not n, `if`(d > max, d-1, d)
end proc
:

@Joe Riel Yes, I agree that it's quite difficult to discover from Maple's help the existence of SetOf, which is the operator that allows a formal Maple property (likely one with an infinite number of members) to be treated syntactically as if it were a set.

It can also be done without SetOf, but at the cost of requiring a bound variable (something that I try to avoid when possible):

is(x::RealRange(-5, Open(infinity))) assuming RealRange(-5,2);
                              true

is(x::RealRange(Open(-5), Open(infinity))) assuming RealRange(-5,2);
                             false

@Tamour_Zubair You've switched the signs of some of the coefficients in your ttt1 calculation. It should begin

ttt1:= abs((25*k1/216 + 1408*k3/2565 + 2197*k4/4104 - k5/5) - (...));

where the ... represents the part that you had correct already.

@Reshu Gupta The highest-order derivatives occuring in your system for each of the 5 dependent variables are

diff(H(x),x,x), diff(F(x),x,x), diff(G(x),x,x), diff(theta(x),x,x), diff(P(x),x).

The first step towards numerically solving the system is to algebraically solve the system of ODEs for those 5 derivatives as if they were simple variables. That means that every ODE has to have at least one of those 5 derivatives. But your EQ1 does not.