Only considering upto the first loop, you have 3 syntax errors:

1. In Cond, the last item is spelled Thata instead of Theta.
2. D notation rather than diff must be used in initial and boundary conditions. So, the last item in Cond should be 
   D[1](Theta)(1,tau) = 0. The [1] here indicates differentiation with respect to the 1​​​^st variable, which is X.
3. In the loop, the eval command should be
   eval([OdeSys, Cond], Nr= NrVals[j])[].

In the next loop (the generation of the numeric points), the eval needs to be replaced with
Ans[k]:-value, and you need to specify values of X and tau to use. Remove the (j); you can't do evaluations of pdsolve(..., numeric) solutions as if they were ordinary functions. The value when used in this way is documented on help page ?pdsolve,numeric.

Yes, as you said, sometimes there's no feasible solution. In those cases, the solve returns [] (the empty list), which cannot be indexed by [1]. Your procedure needs to be prepared to handle that. I recommend that you change your command with eval and solve to this, where the "..."s represent what you already have:

...:= eval(..., solve(...)[])[];

Now the procedure will always return 9 values, and in the infeasible cases those values will be symbolic.

That's curious. I don't know what causes the problem, nor have I investigated it, because there's an easy workaround by which you can add "global" options to any plot or group of plots. By "global", I mean options that apply to the overall plot rather than the mathematical options that specify how to calculate points and colors for the plot. The command that does this workaround is plots:-display:

display(
    fieldplot3d(
        [0, 0, -y], x= -2..2, y= -2..2, z= -2..2, arrows= `3-D`, fieldstrength= maximal(0.5),
        font= [Times, bold, 16], grid= [4, 4, 4], labels= ['x', 'y', 'z']
    ),  #end of fieldplot command
    labelfont= [Times, bold, 40]  #additional global option
)

The axes labels being "cut off" immediately above is just caused by the transfer of the plot from my worksheet to MaplePrimes. They appear normal in the worksheet.

I could've also considered the options font and labels to be "additional global options", but I didn't because they were working as is.

The display command can also be used to combine multiple plots into one presentation. See help page ?plots,display.

The ".mtx" (aka MatrixMarket) file format is a matrix format, not necessarily a graph format. Since the entire mathematical structure of a graph is encoded in its weight matrix (or adjacency matrix if you don't care about weights), this doesn't present any serious difficulties; you just need the tiny additional step of creating a Matrix when going from Graph to ".mtx" file or vice versa.

The example file that you linked is a sparse representation likely suitable for the vast majority of graphs. It's not necessary that you understand the representaion in order to use this. Only nonzeros are represented. The 1st column is the row numbers, the 2nd is the column numbers, and the 3rd is the weights. As you said, the 3rd column could be all 1s, making it an adjacency matrix.

To go from a Maple Graph G to a ".mtx' file, use

ExportMatrix("my_graph.mtx", GraphTheory:-AdjacencyMatrix(G));

To go in the other direction, use

G:= GraphTheory:-Graph(ImportMatrix("my_graph.mtx"));

The GraphTheory:-Graph command doesn't care whether you pass it an adjacency matrix or a weight matrix; either case will be handled correctly without you needing to specify which one it is.

You have three errors:

1. The syntax of pdsolve(..., numeric) requires you to change {ICs}, BCs to 
   {ICs, BCs[]}
2. You have the expression exp(-1000*(x-l)^2) several times. What is x? I guess it should be y.
3. You have 8 BCs; it wants 6. I removed the 1st and 3rd (arbitrarily).

Making these changes, the pdsolve command itself runs without error. You may still encounter numeric errors when you try to get output. Indeed, I think that it's likely that you will.

I'd use grid= [33, 129]. So, that's 33 r values and 129 theta values. It takes 50 seconds on my computer:

CodeTools:-Usage(
    plot3d(
        [A,B,C], r= 1/5..4/5, theta= -Pi..Pi, view= [(-8..8)$3],
        shading= zhue, grid= [33, 129], style= surface, axes= none
    )
);

You have Vectors named err1 and err2 to hold your errors. The command LinearAlgebra:-Norm is akin to an "absolute value" for vectors---it's a way to assess their magnitude as a nonnegtive real number. So, use 

if LinearAlgebra:-Norm(err1, 2) < tol then ...

The 2 indicates the 2-norm (aka Euclidean norm). Other common norms are 1 (sum of abs of elements) and infinity (max of abs). See ?LinearAlgebra,Norm.

Like this:

local D:= Matrix(A, shape= diagonal);
(L,U):= map2(Matrix, A-D, shape=~ [triangular[lower], triangular[upper]])[];

By setting

Digits:= 50;

and adding option complex to your fsolve command (and no other changes), I get this solution in under 2 seconds:

{A = -2.7553365135418814642586082436429575890825402826031,
B = -0.70285804987973303586180028708027467941012949957141}

Note that this solution is real and didn't require initial specification of intervals.

Using the 3rd-party package DirectSearch (available for free download from the Maple Applications Center), and using its command SolveEquations with option AllSolutions, I get 2 solutions: the same one that I show above, and the one reported by Rouben. I'm not confident that this program, or any numeric method, can guarantee the number of solutions to a system of more than 1 non-polynomial equation, even when all univariate restrictions are meromorphic (as is the case here).

@Rouben Rostamian  The procedure can be shortened like this:

tangent_plane:= P-> [P, (sign*primpart@sort@`.`)((<1,3,5>, <x,y,z>)-~'<P>') = 0]:
tangent_plane~(L);

@Laurenso: For your purposes here, the single command sign is equivalent to your signum@lcoeff, but is faster and more robust. It is designed for this type of operation (exact work on polynomials); signum is not.

Here is my version, which uses lists of integers rather than arrays. It's a very short code. The run time is just slightly longer than dharr's.

restart
:
nextone:= L-> local i:= 2, j, n:= nops(L);
    ["", (do for i from (j:= i) + 1 to n while L[i]=L[j] do od; i-j, L[j] until i > n)]
:
LookAndSay:= (n::posint, start::nonnegint)-> local L:= ["", start];
    <start, (to n-1 do parse(cat((L:= nextone(L))[])) od)>
:
CodeTools:-Usage(LookAndSay(50,1));
memory used=309.93MiB, alloc change=27.45MiB,
cpu time=3.52s, real time=3.52s, gc time=281.25ms

Regarding the timing of your Maple code (the StringTools-based code in your Question): You can't blame the compiled StringTools code. The vast majority of the time is spent doing garbage collection at the Maple kernel level, not in the compiled code. This is due to the large number of temporary substrings that you create.

You just needed some minor syntax corrections:

R:= plots:-implicitplot(
    [x = -3/4, x = -1/3, x = 1/3, x = 3/4, y = -3/4, y = -1/3, y = 1/3, y = 3/4], 
    x = -3 .. 3, y = -3 .. 3, color = red, scaling = constrained, gridrefine = 3
);
z:= x + y*I:
plottools:-transform(unapply(evalc([(Re, Im)(1/z)]), x, y))(R);

The circles don't show as complete circles because the lines are not complete (infinitely long) lines; they're line segments.

To get a single solution, you just need to make 3 trivial changes:

1. Correct beta_1 = 0.2 to beta_1:= 0.2;
2. Correct all pi to Pi;
3. Change solve to fsolve.

If you need all solutions, or all real solutions, it's a little trickier.

The example code that you show

f[N1_]=sum[1/(n^3*sin^2[n]), {n,1,N1}];
DiscretePlot[f[x], {x,0,400},PlotRange->All];

is Mathematica, not Maple.

For Maple, you can do a continuous plot of it as

plot(eval(sum(exp(-lambda)*lambda^n/n!, n= 0..N1), lambda= 15.4), N1= 0..50);

If you insist on a discrete plot, simply replace N1 with trunc(N1) inside the sum:

plot(eval(sum(exp(-lambda)*lambda^n/n!, n= 0..trunc(N1)), lambda= 15.4), N1= 0..50);

You need to set your working directory to something that your OS gives you permission to write a file to. This can be done with the currentdir command. For example, on my computer, this works:

currentdir("/users/carlj/desktop");
ExportMatrix(matlabData, A, target= MATLAB, format= rectangular, mode= ascii);