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

@Muhammad Usman Do you agree or disagree that there is something wrong with your value of epsilon​​​​​? 

@Fereydoon_Shekofte Coding Theory is the branch of mathematics that studies error-detecting and error-correcting "codes", which are means of encoding digital/binary information with extra bits (in some optimal way) so that it's possible to detect most accidental^[*1] errors arising from the transmission or storage of the data via potentially noisy media. With some more extra bits, it's possible to correct some errors also.

A prominent example is compact disks (CDs). The code that they use can perfectly correct up to 15%^[*2] of corruption of the disk. That's why dirty and scratched CDs are often still playable, producing exactly the same sound as when they were new.

[*1]  I emphasized "accidental" because, AFAIK, Coding Theory makes no attempt to thwart intelligent, malicious hackers who know the code from getting corrupt data through a channel by design.

[*2]  Getting the 15% maximum correction depends on the corruption being evenly spread over the CD. This is why if you wipe a CD with a cloth to clean it, you should do it radially. A scratch along an arc concentric to the disk can corrupt a large contiguous swath of data and not be correctable. Even if it's much less than 15% of the total disk, it may still be more than 15% of the contiguous data blocks that have been damaged.

What happens if you go to menu File => Restore Backup? If backups are enabled, that should take you to a standard file-selection dialog.

@mmcdara The only place where exp functions and the numeric literal 1000 occur are in the ICs. This is precisely what you did change.

(In the following discussion, I've changed the OP's lowercase l to L because the l is easily mistaken.) 

You seem to be trying to compensate for the OP's very wrong huge epsilon value (> 10^7) by making the argument to exp of much smaller magnitude. IMO, that can't possibly work in any meaningful way----and here's why:  We're given that y is bounded by 0 <= y <= 1, and L is some fixed value of y  (L = 1/5 is used in the OP's example). The piecewise expression that uses epsilon in the ICs is of the (logical) form

• "If y is within epsilon of L, then return a very small positive number, otherwise return 1 minus the same very small number."

So how can an epsilon > 1 be meaningful? I'm not saying that it'll be erroneous in any arithmetic sense; rather, I'm saying that it can't serve any practical purpose. That value of epsilon must be the result of an error the OP made in setting up the problem. 

I do think that the -1000 used by the OP will be problematic, as I mentioned earlier, and its magnitude should be reduced. But I wouldn't go so far as to change it to -1/1000 as you did.

@Muhammad Usman Your value of epsilon is 3.12 * 10^7. That's way too big. I expected epsilon to be something like 10^(-3) or less.

Also, note that in floating-point computation 1 - exp(-1000*(y-l)^2) is often likely to be identically 1 (a phenomenon known as machine epsilon), which (I'm just guessing here) is precisely what you're trying to avoid.

@acer I think that your disparagement in the ending paragraph of your Answer may be due to conflation of the name-forming quotes `...` and the empty-named function ``(...). Only the former can lead to inappropriate changes in typography, such as italics. Here is each applied to the OP's situation (assuming one wanted to directly input the expression as inert). Compare the 2D output of these:

1. `3`*(a - 1)  #3 becomes an italicized name.
2. 3*``(a - 1)  #produces desired output; a-1 remains type `+`.

The expand command will turn #2 back into what one would've had if `` weren't used in its input.

I think that the OP may want to deconstruct to an inert form an already-automatically-simplified expression. In that case she can use (assuming expr is an expanded polynomial such as 3*a - 3):

content(expr) * ``(primpart(expr))
or (same thing in a more-functional style)
(content * (``@primpart))(expr) 

Vote up.

I converted this from a Question to a Post yesterday. I only mention it so that it'll appear on the "Active Conversations" list. The conversion process seems to remove it from that list.

The default grid is [49,49], so your [200,200] is more than 16 times larger, so it likely takes more than 16 times the amount of time.

@NIMA112 In  your case, singularities occur along the lines in the A-B-plane where the arguments to your lns are 0. Specifically, those lines are B = 0, A = -B, and A = -100*B.