@Adam Ledger Could you please write those formulae in a more-human-readable format rather than Latex or MathJax? I know that that incurs a risk of a loss of precision. That's a risk that I'm willing to accept.

By "first term of the asymptotic expansion of f(n)" what's meant is some relatively simple function a1(n), preferably eventually^[1] monotonic and analytic, such that limit(f(n)/a1(n), n= infinity) = 1. Note that the PNT does nothing more than say that the first term of the asymptotic expansion of the prime-counting function pi(n) is n/ln(n).^[2] So, just getting one term can be quite an achievement.

[1] "Eventually" means for all sufficiently large n.

[2] The function Li(n) = Int(1/ln(x), x= 2..n) can also be used for the first term, and it gives a more accurate approximation.

@Adam Ledger I am not an expert on the Prime Number Theorem, but perhaps I can give some help anyway. I read your StackExchange post that you linked to. 

So, you have a function defined on the positive integers, which I'm going to call f(n) (feel free to change the f). Here it is:

q:= n-> min(divisors(n*floor(ithprime(n)/n) minus {1}));
f:= n-> (n*igcd(n, n*q(n))-floor(sqrt(n))*igcd(floor(sqrt(n)), q(n)))^(1/2);

Please confirm that I have that right.

So, you want an asymptotic approximation to that. The evidence from the StackExchange post's replies suggests that the first term of that will be n.

@Gabriel samaila For the sake of completeness and full disclosure, here's the code that does the check that I referred to in my previous Reply:

_V0:= <V0, 0.5, 0.75, 1, 1>:
_Pr:= <Pr, 0.7, 0.7, 0.7, 0.015>:
_Pm:= <Pm, 0.5, 0.5, 1.0, 0.5>:
_Ha:= <Ha, 5, 5, 5, 0.5>:
T1:= <_V0 | _Pr | _Pm | _Ha >:
_Cf:= Vector(5, [C__f]):
#######
#This little hack is required to set Nb to 0 because it's a 
#denominator in ODEs. By setting Nt to 0 first, the fraction with Nb in the
#denominator disappears:
savODEs:= copy(ODEs):
ODEs:= eval(ODEs, [Nt, Br, Sc]=~ 0):
#End of hack.
######
for k from 2 to 5 do 
   sol:= Solve(convert(T1[1,..] =~ T1[k, ..], list)[], ([Nt, Br, Nc, Nb]=~ 0)[]);
   _Cf[k]:= eval(C__f, sol(.5)) #.5 is arbitrary between boundary points.
od:
#Restore ODEs to its original pre-hacked state:
ODEs:= savODEs:
#######
<T1 | _Cf>;

I converted your Post into a Question. I also needed to lengthen the title because there's a 15-character minimum.

@Carl Love Okay, here is a fair test that shows that iterating with nextprime is many times faster than iterating with ithprime when the primes are in the range just above the limit of Maple's fresh-out-of-the-box pre-stored knowledge of prime ranks, which is where you're currently getting stuck. To generate the first 30,000 such primes, nextprime is 25 times faster. But, the slowness of ithprime will increase (quadratically!) as you proceed due to an idiosyncracy of the way that Maple handles lists.^[1]

Download nextprime_vs_ithprime.mw

[1] Here are some details of that idiosyncracy: The ordinates (or indices) of the remember table of procedure `ithprime/large` are stored in the list `ithprime/global/list`. Use showstat(`ithprime/large`) to display the code of the procedure. On lines 6, 10, 17, and 22, you can see that items are added to the list as it acquires new knowledge about the ranks of larger primes. It is well known that adding items to long lists in Maple is extremely inefficient; indeed, the time required is proportional to the length of the list; so, using ithprime in a loop requires time at least proportional to the square of the number of iterations.

You can also see on lines 3-4 of the procedure that ithprime simply iterates with nextprime anyway when it's operating beyond the limit of its stored knowledge. That alone proves that nextprime must be at least a little faster for those primes. But the list idiosyncracy makes nextprime much faster.

@David Sycamore Rather than make a comparison of ithprime and nextprime (because I don't have time at the moment to construct a fair comparison of the two), I made the following test, which shows that for the range of primes that you're working with (I think), the nextprime method iterates through them at more than 30,000 primes per second  on my machine:

restart:
P:= proc(P::posint, n::posint)
local p:= P;
   to n do p:= nextprime(p) od;
end proc:
CodeTools:-Usage(P(3*10^8, 3*10^4)):
memory used=173.43MiB, alloc change=0 bytes, cpu time=906.00ms, real time=902.00ms, gc time=46.88ms

@Gabriel samaila If you make Nt=0 first, that will eliminate the term that has Nb in the denominator. After that, you can make Nb=0. Since the four unknown parameters now have numeric values, there's no need use fsolve; you can just use the procedure Solve that I already posted. I've done this, and I still get values for C__f completely different from those shown in Table 1.

@Gabriel samaila I haven't tried it yet. I think that you mean making Nb a small value, not Nt? I will try it.

@tomleslie My experience is that not putting restart in its own execution group is more likely to cause a problem with interface commands than with others. I think that it has something to do with synchronization between the GUI and the kernel.

@Rouben Rostamian 

Good workaround. Vote up.

My guess about what's causing the error (just a hunch) is that it happens when the size of the window needs to be expanded for some reason that Explore didn't originally anticipate. In this case, the superscripted exponents cause a slight expansion in the vertical extent of the plot.

Explore is by far the single largest procedure that I've seen written in Maple. Listing its code into a string with sprintf, the string is 56,591 bytes. It has several magic formulas for estimating the height and width of the window needed.

@samira moradi For example, the following command will find all the roots in the interval 0..99 for all integer B in 1..30. This command takes about 5 minutes to execute:

Lambda:= table(
   [seq(
      j= 
      Student:-Calculus1:-Roots(
         eval(F[lambda[i]], [y= 0.9, k= 0.003, B= j]), 
         0..99
      ),
      j= 1..30
   )]
);

Then Lambda[25], for example, will be a list of the 55 roots for B=25, and Lambda[25][30] will be the 30th root for B=25.

@samira moradi The command Student:-Calculus1:-Roots will give you all the roots within a bounded interval.

@samira moradi Given that mostly trig functions are used, I strongly suspect that there are an infinite number of roots.

@nm You wrote:

• Not that I even know what it means being a newbie in Maple.

As I've told you before, setting kernelopts(opaquemodules= false) makes the syntax A:-B apply to both exports and non-exported locals of module A; whereas the default setting, true, means that it applies only to exports^[1]. I guess that the only significant thing about it that I didn't say before is that that is the only purpose of opaquemodules.

[1] The exception is if A is literally the keyword thismodule. This can only be used inside a module, and then thismodule:-B can refer to both exports and non-exported locals of the current module regardless of the setting of opaquemodules.

@tomleslie Thanks, I didn't know that. That's very useful. I hope that that setting can be saved for future sessions.