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

@SHIVAS I knew that there were more errors; I simply couldn't correct them all at once. To proceed, I need to know the numeric values of X and tau at which to evaluate diff(Theta(X, tau), X) for your desired Matrix output.

@RezaZanjirani The problem with the method that you just proposed is that using it requires you knowing beforehand the values of alpha that are feasible.

Your PDE system has 1 dependent and 2 independent variables. You should correct the title of this Question.

@MACi You seem to be apologizing for Posting, but I don't think that you have anything to apologize for. Your contributions are welcomed! I think your suggestion for a Coding Theory package is a good one. If you have questions while you're learning Maple, please send them. Of course, they should go in the Questions area rather than Posts.

And please don't ever delete appropriate content that has been responded to. It can be extremely offensive to the person who wrote the response. If you want to retract something, post a Reply saying so.

@Fereydoon_Shekofte Thank you; that was so nice of you to say.