@Oak Staff 

Before I explain it, I need to make a correction. There was a bug in my code that would only turn up when the last die had the same value on multiple sides (a rare situation). Here's the corrected code:

Tallys:= proc(D::list)
option remember;
     table(sparse, Statistics:-Tally(D))
end proc:

DiceSumRec:= proc(S,D)
option remember;
local k;
     if nops(D)=1 then  Tallys(D[])[S]
     else  add(thisproc(S-k, D[2..])*Tallys(D[1])[k], k= indices(Tallys(D[1]), 'nolist'))
     end if
end proc:

DiceSum:= (S, D::list(list))-> DiceSumRec(S,D) / `*`(nops~(D)[]):

Here's some explanation. The denominator at the end of DiceSum is `*`(nops~(D)[]). This is the product of the number of sides of each die. That's the total possible number of rolls. The numerator is the number of those rolls that have sum S. These rolls are counted by recursively examining each roll and counting it if its sum is S.

To explain the recursion in more detail, let's look at the example that I already used: getting a sum of 10 with three standard six-sided dice. We can get a 10 by getting a 1 on the first and a sum of 9 on the remaining two. That can happen 4 ways. Or we can get a 2 on the first and 8 on the remaining two. That can happen 5 ways. Likewise, 3+7 can happen 6 ways; 4+6, 5 ways; 5+5, 4 ways; and 6+4, 3 ways. The total ways is 4+5+6+5+4+3 = 27. The probability is 27/6^3 = 1/8.

If you want more explanation about the nuances of the Maple syntax, let me know.

To MaplePrime Adminstrators: Notice that this Question has a Vote Up (given by me), yet the OP still has reputation 0. This should be impossible.

@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.