If this is an attempt by you to find the longest path, it won't work.

@vs140580 From the Wikipedia article "Hypercube graph", I learned that for distinct vertices u and v in HyperCubeGraph(n) there is a Hamiltonian path between them iff their Hamming distance is odd. Thus, if the distance is odd, then the length of the longest path is 2^n - 1. If the Hamming distance is even, I don't know, but I suspect there's an equally simple formula.

@vs140580 If you'd want to restrict attention to certain highly symmetric special cases, such as hypercubes, then a fast answer is possible. The vertices of HyperCubeGraph(n) can be represented as n-lists of 0s and 1s. Because of the symmetry, the relationship between any two vertices is (up to isomorphism) entirely determined by the number of positions that are different in their 0-1 representations (known as their Hamming distance). Thus, there are only n possible relationships between pairs of distinct vertices. Thus there are only n possible values for the length of the longest path.

@sand15 The issue is that Sample uses an hfloat array by default. If you map those integerizing functions over an hfloat array, the result is still an hfloat array, and thus can only contain hfloats. This allows those commands to work under evalhf. But if you use an sfloat array, as you did, then mapping those functions changes the datatype of the array, because now there's no worry about evalhf.

But if you first convert to a list, as I did, before mapping the integerizing function, it works fine. As my code shows, this is many times faster than using an sfloat array.

You asked:

• Couldn't  it be interesting that these three functions also operate on hfloats?

Yes, it is, and exploiting that fact can lead to remarkable speed gains (20X - 30X). Hardware floats can exactly represent and compute with integers under 2^50 (or so). This is one reason why the LinearAlgebra:-Modular package is so fast. For the functions that can work with evalhf, see ?evalhf,fcnlist.

@vs140580 The longest-path problem for general graphs and digraphs is known to be NP-hard, meaning that there's provably no fast algorithm for it. You need to check every path.

@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));