@mmcdara The Algol that you quoted

i:= k; for i: = i— 1while co x eps2 < p A i > 0 do

translates to Maple as

for i from k-1 by -1 to 1 while co*eps < p do

Although I haven't read through it super carefully, I don't see any typos in the Algol code.

@Mo_Jalal It would've been a worthy experiment with respect to the run time if you had restricted the computation to -20..20. But that's not what your code shows. It shows that you computed over -80..80 and then plotted over -20..20. The time to plot the points once they're computed is nearly zero.

Gonzalo: It will be trivial, although tedious if you're not a good typist, for YOU to type the three pages of Algol code given in the linked PDF into a plaintext file or a Code Edit Region in a Maple worksheet. After that, it should take an hour or less for one of us to convert it to Maple. The code is quite straightforward, and translating it to Maple should be trivial.

@Mo_Jalal The plot of any continuous function can be done by precomputing some points, like you've done. I don't see any benefit to doing it that way. 

You've computed points on -80..80 but only plotted on -20..20. Is there a reason for that discrepancy?

@Reshu Gupta Change the dsolve command to this:

sol:= dsolve(EQ union IC, numeric, maxmesh= 2^16, initmesh=512, abserr= 0.4e-3, method= bvp[trapdefer]):

The plots will take significantly longer to produce (1-2 minutes), but you will get them.

@Carl Love And here's an idea for representing the terms more compactly.
 

Download PolyDeconSorted.mw

@ogunmiloro In the initialization you defined C[f], but in sse you used C__f. Of course, both usages must be spelled exactly the same.

You've ignored advice that I gave you for an earlier Question about this same problem: Use option maxfun= -1 in the dsolve command. This will avoid all the "cannot evaluate" errors.

Making these changes, the optimization command will return with "no improved point could  be found." I can't help you with that at the moment. It's not due to a syntax error. Rather the problem either requires some fine tuning or it's intractable.

@Bohdan Here it is:
 

Download PolyDeconSorted.mw

@Bohdan I think that you may have misunderstood my question about the rhos, so I'll try to clarify it. I know that all the rhos will appear in every numerator. My question is If the numerators were expanded, would there be any **single term (aka a monomial)** with two or more rhos multiplied together?

Regarding numer: First use numer, then collect, then coeff. I'll give a detailed example in a later Reply.

Regarding the display of the sum that you want: This should be handled by Physics:-Library:-Add, but I can't get it to work in this case (see ?Physics,Library). Someone who knows the Physics package better may be able to help. However, if you want the computation instead of the display, I can easily help you with that.

Regarding sorting: I'll give a detailed example in a later Reply.

@Mo_Jalal 

Okay, if the upper limit is changed from infinity to a real constant, then the lower being nonreal isn't a problem.

The sum is a problem. For reasons that I don't know, I get an error "singularity encountered in the interval of summation". I can't see any singularity. Anyway, I found a way to work around this problem, shown below.

By using numeric integration with an accuracy of 5 digits (which should be sufficient for any plot), I can get a 50-point plot in 5 minutes. In the code below, I've highlighted my changes in green:

restart:

assume(d,real):
b:=-gamma1*tau0+I*tau0*Delta-2*a*(k-1)*xr:
a:=1+I*c;
c1:=sqrt(conjugate(a))/tau0:
c2:=0.5*((conjugate(b)/sqrt(conjugate(a)))):
lambda1:=2*(1-I*d)*(t-k1)/j+(1-I*d)^2:
lambda2:=t*(1-j)*k1/j+1/sqrt(2)*(1-I*d):
lambda3:=c1*(t-k1)/j+c1*(1-I*d)+c2:
J1:=sqrt(Pi)/sqrt(2)*(1-erf(lambda2)):
J1mod:= abs(J1)^2:
g1:=0.5*sqrt(Pi)*tau0*exp(c2^2)*exp(-conjugate(a)*((k-1)*xr)^2)/(sqrt(conjugate(a))):                       
F2:=-sqrt(2)*Int(exp(-x^2)*erf(sqrt(2)*c1*x+lambda3),x=-lambda2..infty, 'epsilon'= 1e-5);
J2:= Sum(g1*J1+g1*F2, k= 1..K);
J2mod:= abs(J2)^2:
W:= unapply(no*J1mod+Omega^2*J2mod, d):
Omega:=0.01:no:=1:Delta:=0:tau0:=1:c:=0:gamma1:=0:j:=1:k1:=0:t:=2*Pi:xr:=1:
infty:= 100: K:= 1:
W:= subs(Sum= add, eval(W));
CodeTools:-Usage(plot(W,-50..50,numpoints=50,axes=boxed,title=tit,color=black,font=[2,3,18],thickness=2,tickmarks=[3,2],titlefont=[SYMBOL,14],font=[1,1,18],linestyle=1));

@Arif Ullah khan I sent you an email about 24 hours ago, but I haven't received a reply.

@Reshu Gupta I didn't mean to discourage you. You're welcome to post that boundary-layer Question as a new Question.

@Carl Love As can be seen from my code above (which you've confirmed in another Reply is correct), your entire 12x12 system can be generated from 7 fairly easy and quick-running lines of Maple code. With one more line, I can extract the coefficient matrix and right-side vector of the system:

(A,B):= LinearAlgebra:-GenerateMatrix(sys, Ps);

With some simple adjustments, which wouldn't increase the complexity of the code in the slightest, it could be used to generate the 27x27 matrix that you initially inquired about. So, all the effort and back-and-forth dialog about importing the matrix could've been replaced by a few simple lines of Maple code. (I'm not saying that you should've known this at the time.)

The take-away message is: Consider the possibility of a pure-Maple solution first.

@Bohdan I have a very low personal psychological tolerance for GUI sluggishness, and your worksheet "Problem_3.mw" for the12x12 system slightly exceeds that tolerance. Would it be possible for us to "practice" similar simplification operations with a smaller system, say the 8x8 system that you referred to earlier (if it's full rank)? There are easy and very fast commands that can replace all the operations that you laboriously do with the mouse. One such command is numer, which in this case will extract the numerator polynomial of a rational function. (I think that you realize already that all the denominator polynomials are the same, being the determinant of the coefficient matrix of the original system.) I'd be happy to teach you these commands on a smaller system.

Is this statement correct: No term in any of the expanded numerators contains more than one rho?

My guess is that you want to express these numerators in the form 

Sum(rho[i,r]*Sum(Product(pi[...], ...)*Product(gamma[...], ...), ...)

Is that correct?

The hard part is seeing the patterns of the indices. I don't think that simplify will help you see that. What might help, though, is having the terms sorted by those indices. Do you agree? This can be easily done.

There is a bug in Maple's reporting of the memory allocation. The amount that the operating system says is allocated is almost always (perhaps always?) exactly double the amount that is reported by those "alloc=..." messages in the command-line interface. In the GUI interfaces, that same information is reported as "Memory: ..." on the bottom line of the window. I have noticed this bug for years.

So, the message "alloc=32540.2MB" indicates that the program was using about 64GB (your machine's limit) at the time of the crash.