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

Your Question is a very common situation with this kind of numeric BVP, known as a boundary-layer problem. It's very important to pay attention to the distinction between these two error messages, which only differ by one word:

1. Newton iteration is not converging;
2. initial Newton iteration is not converging.

Several solving tips:

1. Unlike a lot of error messages, these two do not indicate that the user has entered something wrongly. Rather, they indicate features inherent to the problem itself. And they're not just mathematical features: "Boundary layers" are things that really exist in the moving fluids that these BVPs model.
2. These can be very difficult problems. They are one of the most-common sources of Questions on MaplePrimes. Unfortunately, almost no-one here calls them "boundary layer", so that term doesn't help with searches. But the vast majority of Questions about numeric BVPs are about this kind.
3. Read, and reread, the help page ?dsolve,numeric_bvp,advanced. Even an expert with decades of experience such as myself rereads some help pages many, many times over the years. 
4. The continuation option can only help with the "initial" case; indeed, it's the only purpose of the option, and, conversely, the option is usually the only thing that can help.
5. Because of the existence of that option, the "initial" case is relatively easy to fix compared to the non-initial case.
6. Fixing the "initial" case often leads immediately to the non-initial case. But using the continuation option was not in vain because the "initial" problem needs to be fixed before you can try to fix the non-initial one.

Nearly all of these BVPs have a large number of parameters. Five years ago, I wrote an extensive worksheet and Post about this, "Numerically solving BVPs that have many parameters". Several of Maple's numeric IVP solvers allow for symbolic parameters whose values do not need to be specified at the time dsolve is called. The BVP solvers do not. The problem in that worksheet is a boundary-layer problem, but the purpose of the Post and worksheet is to show an easier way to deal with parameterized BVPs. It's purpose is not to show how to deal with these "Newton convergence" issues. Using the techniques in that worksheet, you can produce plots where the independent variable(s) are the parameters.

After correcting the "initial Newton..." errors, there are usually still combinations of parameters that are physically meaningful but lead to the non-initial "Newton" error. This is often because there are multiple solutions for those parameter values^[*1]. For this you need to use the approxsoln option, also described on ?dsolve,numeric_bvp,advanced. The approximate solution used can be that from a nearby group of parameter settings for which convergence was obtained.

[*1] Unlike IVPs, solution uniqueness is usually not guaranteed for BVPs. A simple example is

{diff(y(x), x$2) = -y(x), y(0)=0, y(Pi)=0}

which has the solutions

{y(x) = 0, y(x) = sin(x)} 

and others.

This isn't related to your primary Question; I just want to correct a nuisance from my original Answer. That cutting off of the axes labels also showed in the worksheet, although it wasn't as severe as shown in the Answer, causing me to not notice it. It can be prevented by including 

orientation= [45, 75, 3]

in the options passed to display. The numbers are axial rotations in degrees.