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

@Tamour_Zubair According to the Wikipedia article "Runge-Kutta-Fehlberg method", the formula for h_new is

h_new:= 0.9*h*(tol/ttt1)^(1/5);

I don't think that your using 0.84 instead of 0.9 makes much difference, but I think that having h inside the parentheses does.

You have 5 dependent functions (F, G, H, P, theta) and 4 equations (EQ2, ..., EQ5). Where is EQ1?

You also have a trivial syntax error: You've attached "= 0" twice to each equation. But the important issue is the missing equation. 

@dharr There's a Wikipedia article "Pseudospectrum" about it. The pseudospectrum of a square complex matrix A is the set of all eigenvalues of A+E for all matrices E of sufficiently small norm. "Pseudospectrum" is a bad name; I'd call it the "fuzzy spectrum". I think that plotting it involves drawing contours around the eigenvalues of A.