What is the range of xi that you want to use? If you tell me that, I think that I can show you how to plot it.

You have some small syntax problems. You will confuse Maple by trying to use both p and p[0] in the same program. You should change p[0] to p__0. Likewise for q[0] and q[s].

@Preben Alsholm 

So the plot that I was able to get with Re is actually bogus---some noise from round-off error.

@Carl Love 

Here's a complete prime factorization program that calls the Pollard repeatedly. It should work completely on any composite upto 12 digits. If it encounters a composite for which Pollard returns FAIL, that number will appear in square brackets in the factorization.

Pollard:= proc(
     n::And(posint, Not({identical(1), prime})),
     {B::posint:= max(iroot(n,4), 2^20)},
     {A::And(posint, Not(identical(1))):= 2}
)
description "Pollard's p-1 algorithm";
local q:= 1, lgn:= ilog2(n), g:= igcd(A,n), a:= A;
     if g > 1 then  return g  end if;
     while q <= B do
          q:= nextprime(q);
          a:= a &^ (q^trunc(lgn/ilog2(q))) mod n;
          if a=1 then return thisproc(n, ':-A'= A+1) end if;
          g:= igcd(a-1,n);
          if g > 1 then return g end if
     end do;
     FAIL
end proc:

IFactor:= proc(N::And(posint, Not(identical(1))))
local Q:= SimpleQueue(N), PF:= 1, n, f;   
     while not Q:-empty() do
          n:= Q:-dequeue();
          if isprime(n) then  PF:= PF*``(n);  next  end if;
          f:= Pollard(n);
          if f=FAIL then  PF:= PF*``([n]);  next  end if;
          Q:-enqueue(f);
          Q:-enqueue(n/f)
     end do;
     PF
end proc:

@Declan 

Could it be in your assignment that you can assume the n has exactly two prime factors? That would be typical of a cryptography application. You should use Maple's isprime to detect whether any generated factor is prime. If it isn't, then you can run Pollard again on it. I made a small modification to the posted procedure so that it's guaranteed to find a proper factor of any odd composite with fewer than nine digits. You can increase this guarantee by increasing the 10^4.

@Sunny9 Use op(1, L1fh).

This isn't a complete answer, but certainly you want to avoid the jpeg format for anything with crisp lines or text (that includes tick labels). Try GIF.

The plot shows g bounded below by 0, and equal to abs(1+d) for all c greater than a certain value.

Acer and VV, thanks for clearing up my lack of knowledge about the efficiency of @. Acer, your first group of three tests is to me definitive. The second group of three tests is insignificant regarding the efficiency of @.

See ?CurveFitting,ArrayInterpolation.

@vv I'd be very interested in discussing the speed of @. I chose what I did because it's easy to map without creating a user-defined function. Let LL be a very long list of short lists of digits. Which of the following is fastest?

map(parse@cat@op, LL);

map(`@`(parse,cat,op), LL);

map(L-> parse(cat(L[])), LL);

I'm prejudiced about the speed of the third because it uses a user-defined function. But perhaps my prejudice is unjustified.

Do you have some method to do it "by hand"? Perhaps by analyzing that we can figure out a way to do it with Maple. Of course, any method that starts with the final form is invalid. Maple has no trouble showing that the final form is equivalent to the starting form.

@mskalsi Yes, I did read Tom's solution. It's a fine solution. I simply propose an alternative.

@JessyOw My point is that whatever you want to add to each term, you can do it with +~.

I've told you before that Psi isn't a constant. You're taking code that I wrote regarding gamma and misappropriating it to Psi.

@Markiyan Hirnyk In your code for the first plot in this Answer, you plotted only 53 of the 182 points returned by DirectSearch. Why? Was it that the residuals for the other points were too high?