@Christian Wolinski

collect(evala(Norm((lhs-rhs)(A))),Y,`convert/parfrac`);

@Carl Love 
Yes, I have corrected.

@Carl Love 

No, I din not think much about it. But it seems that a more general fact is true.
If p is a fixed prime and 2 <= k <= p-2 is also fixed then the events
A = {a in {1,...,p-1} :  a is prime} and
B = {a in {1,...,p-1} : a^k mod p  is prime}
are "asymptotically independent". For k=p-2 one obtains the inverses mod p.

@Carl Love 

Actually the two events are not independent (for a fixed p). They seem to have some kind of "asymptotic independence" (p --> oo).
It's interesting your intuition about this!

@acer 

It could be that the integral is exactly 1/8 but I have not time now to check.

@digerdiga 

The numerical integration is easier for rectangular domains and it is always possible to reduce the integral to such domains via a change of variables.
Note that we may include the domain D in a larger rectangle R and define the function 0 in R \ D but this is not the best solution because the new function is (generally) not continuous.

@brian bovril 

This version of alphametics is much faster and uses an intelligent algorithm (not the brute force one). It works even for systems of "equations".
The author is Robert Israel.

@Gabriel samaila 

I see that you have chosen infinity = 1. If you choose a larger value the results are closer.
But I am a mathematician, not an engineer; you will have to fine tune it yourself; for a too large value dsolve may fail or must also be fine-tuned (there are many options for this).

@baharm31 

It works for me. After all it's a linear ode system with constant coefficients which is well known to Maple.

2DOFnumerical-x.mw

@Gabriel samaila 

You must also insert a * (or a space) after M as Rouben noticed. Note that in Maple  7(x) simplifies to 7, not to 7*x.

@brian bovril 

1. Have you tried to use the icon "Open the current help page in a worksheet window" and run from there?

2. Does the Maple compiler work on your system?

@ola123 

Actually it's easier to solve the system by hand. For a,b,c,d complex the solutions are:

{a=0,b=0}, {b=0,c=0}, {b=0,d=0}, {a=0,d=0}, {c=0,d=0}, {a=0,c=0}, {a=k*c, b=k*d, where c,d in C\{0}, |k|=1}

As you see, even forcing a>0 there are solutions with b,c,d complex (nonzero imaginary parts).

@Joe Riel 

I like the program as it is and I'd leave the modifications for an interested user as exercises.
Another exercise could be the reverse problem: starting with a numeric equation such as 29786+850+850=31486 , find a puzzle for it. Of course a large list of words will be needed here; maybe one of the words could be prescribed.

@Christian Wolinski 

Yes but the original polynomial has practically the degree 8 (the K[4] factor is obvious).

@Christian Wolinski 

Why dou you say it's simpler? It has 4 indeterminates, more terms and almost the same degree.
I think that Factor() mod 2 was not used much for >2 indeterminates.