@Carl Love The following is an implementation, essentially, of your originally posted algorithm. Note that there's no need to explicitly compute the factorials because, for example, ((2^2)^3)^4 = 2^(4!).

Pollard2:= proc(
     n::And(posint, Not({identical(1), prime})),
     {B::posint:= max(iroot(n,4), 2^20)},
     {A::And(posint, Not(identical(1))):= 2}
)
local a:= A, g:= igcd(a,n), i;
     if g > 1 then  return g  end if;
     for i from 2 to B do
          a:= a &^ i mod n;
          if a=1 then  return thisproc(n, ':-B'= B, ':-A'= nextprime(A))  fi;
          g:= igcd(a-1,n);
          if g > 1 then  return g  end if
     end do;
     FAIL
end proc:

If you had two more algebraic equations (no derivatives needed), then the system could be solved as a system of DAEs. Can you come up with two algebraic equations describing x(t) and y(t)?

@Declan Could it be that your professor has given you an n which can't be factored by Pollard's p-1 and hopes that you will figure out why it can't be factored by that method? I believe that your n can't be factored by Pollard's p-1 in any reasonable amount of time. Here's my argument:

I used Maple's ifactor command to factor n. This takes less than a minute. As I guessed (for a reason that I gave in an earlier Reply), n has exactly two prime factors, which I'll call p and q. Then I used ifactor again to factor p-1 and q-1. The largest prime factor of q-1 has 27 digits; the largest prime factor of p-1 has 17 digits. Call this latter prime r. The number of primes less than r can be very closely approximated by

evalf(Int(1/ln(x), x= 2..r));

This is the number of iterations of the loop in the algorithm that would be needed to find a factor of n. If we can do one million loop iterations per second, it'll take 25.6 years to do that number of iterations.

Regarding the differences in our algorithms: Both algorithms will reach the same conclusion for any given n. Yours just takes much longer because of many wasted loop iterations (there's no need to check the gcd when i+1 isn't prime). My version of the algorithm is also described in this very-easy-to-read expository paper: http://www.math.mcgill.ca/darmon/courses/05-06/usra/charest.pdf

@lassa The only other important change that I made to your code was to change the variable of integration from q to qs in procedure G. The unimportant changes were that I removed superfluous parentheses and used my preferred spacing.

@marc sancandi All you have to do is resave MC as a text file. Just use a save command, making the filename end .txt instead of .m. Then you can attach that file to a post here. 

No experiment is valid with irreproducible results. Most readers here want to reproduce your results before commenting. It has nothing to do with feeling attacked.

@Carl Love Well, does the above do the job for you? 

@tomleslie 

Maple's numerical PDE solvers are restricted to systems with two independent variables. The symbolic PDE solvers can handle more.

Please post your worksheet as an attached file. Your plaintext code above is too voluminous and messy to copy-paste. 

@John Fredsted Yes, that wasn't my intention.

First, I need to filter the output processing also. So, I change the procedure to

`&B`:= (ex::uneval)->
     subsindets(eval(subsindets(ex, integer, convert, decimal, 2)), integer, convert, binary):

Now you can do

&B `+`(11, 10, x);

You can replace integer with numeric in my code. But I thought that integer would satisfy the OP's purposes.

@John Fredsted 

`-`(2,1,3,4) is the same as 2-1-3-4. `/`(2,1,3,4) is the same as 2/1/3/4. Both infix variations correspond to the way that I would enter such expressions, either in a program or in a hand calculator. I eschew extra parentheses.

Here's a variation that more gracefully handles situations where there are non-integer extra arguments:

`&B`:= (ex::uneval)-> convert(eval(subsindets(ex, integer, convert, decimal, 2)), binary):

(Using evalindets instead of eval@subsindets produces weird, incorrect results. I don't know why.)

@tomleslie 

Yes, that's the post that I meant.

Do you simply mean that you want a procedure to shift the index variable? I believe that I posted a procedure for that a few months ago. Unfortunately, I'm not very good at using the search facility on MaplePrimes. Someone else may find it and post a link.

@Carl Love 

Here's a completely different approach to the dice-sum-probability problem using generating functions. The calling sequence and restrictions are identical to the procedure above.

DiceSumProb:= proc(S, D::list(list))
local e,d,x;
     coeff(combine(expand(mul(add(x^e, e= d), d= D))), x^S) / `*`(nops~(D)[])
end proc:

If you want the entire probability distribution for a particular set of dice, it'd be more efficient to generate it in one fell swoop than to call the above for each value of S. Let me know if that's what you want.

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.