Here is a problem from SEEMOUS 2017 (South Eastern European Mathematical Olympiad for University Students)
which Maple can solve (with a little help).

For k a fixed nonnegative integer, compute:

Sum( binomial(i,k) * ( exp(1) - Sum(1/j!, j=0..i) ), i=k..infinity );

(It is the last one, theoretically the most difficult.)

Sudoku is a well known Latin square type game, see https://en.wikipedia.org/wiki/Sudoku

Here is a Sudoku game and its (unique) solution:

A,Sol:=  # A = Sudoku matrix, 0 for each empty cell
Matrix(9, [
[0,0,3,0,9,0,1,0,0],
[0,5,0,3,0,0,7,0,0],
[1,0,2,0,0,5,0,6,4],
[0,1,0,0,2,0,9,0,0],
[2,0,0,6,0,3,0,0,1],
[0,0,7,0,8,0,0,3,0],
[7,6,0,9,0,0,8,0,5],
[0,0,8,0,0,7,0,9,0],
[0,0,4,0,6,0,2,0,0]]),
Matrix(9, [
[4,7,3,2,9,6,1,5,8],
[8,5,6,3,4,1,7,2,9],
[1,9,2,8,7,5,3,6,4],
[3,1,5,7,2,4,9,8,6],
[2,8,9,6,5,3,4,7,1],
[6,4,7,1,8,9,5,3,2],
[7,6,1,9,3,2,8,4,5],
[5,2,8,4,1,7,6,9,3],
[9,3,4,5,6,8,2,1,7]]);

The procedure which follows is a very compact Sudoku solver. It uses Groebner bases. I hope that you will like it.
The input is the Sudoku matrix and the solution matrix is returned.
Note that the Sudoku matrix must be valid and must have a unique solution.
(Otherwise, theoretically, the error "Invalid Sudoku matrix" should appear.)
Note also that the procedure may be very slow for some games or Maple may crash. This happened to me once with a very "hard" matrix.

I was impressed that Maple's implementation for Groebner bases works now so well for this problem!

A few years ago on this site: http://www.mapleprimes.com/questions/131939-Calculating-Groebner-Basis-For-Sudoku
it was an attempt to solve the problem with this method but it failed (due to wrong number of polynomials).

sudoku:=proc(A::'Matrix'(9,integer))
local x_A,x,Q,R,r, i,j,u,v,G;
Q:=proc(X,Y) normal((mul(X-i,i=1..9)-mul(Y-i,i=1..9))/(X-Y)) end;
x_A:=seq(seq( `if`(A[i,j]>0,x[i,j]-A[i,j],NULL),i=1..9),j=1..9);
R:={seq({seq(x[i,j],j=1..9)},i=1..9), seq({seq(x[i,j],i=1..9)},j=1..9),
    seq(seq({seq(seq(x[3*u+i,3*v+j],i=1..3),j=1..3)},u=0..2),v=0..2)};
G:=Groebner:-Basis({seq(seq(seq(Q(u,v),u=r minus {v}),v=r),r=R),x_A},'_vv');
if nops(G)<>81 then error "Invalid Sudoku matrix" fi;
eval(Matrix(9,symbol=x), `union`(map(u->solve({u}), G)[]));
end:

sudoku(A) < A; # Solving the previous game

# Let's solve another one:
A:=Matrix(9,9,[[0,0,0,4,0,0,0,8,0],[0,5,2,7,0,0,4,0,0],[3,0,0,0,0,0,0,0,0],[5,1,0,8,0,0,0,0,0],[0,0,0,5,0,0,6,7,0],[0,9,0,0,7,0,0,0,3],[2,4,0,0,0,5,0,0,0],[9,0,0,0,0,0,0,3,8],[0,0,0,0,0,0,9,4,0]]):
sudoku(A) < A;

Matrix   # A Sudoku matrix which crashes Maple!
(9,[[8,0,0,0,0,0,0,0,0],[0,0,3,6,0,0,0,0,0],[0,7,0,0,9,0,2,0,0],[0,5,0,0,0,7,0,0,0],[0,0,0,0,4,5,7,0,0],[0,0,0,1,0,0,0,3,0],[0,0,1,0,0,0,0,6,8],[0,0,8,5,0,0,0,1,0],[0,9,0,0,0,0,4,0,0]]):

In this post I want to present an easy method to obtain a discrete parametrization of a surface S defined implicitly (f(x,y,z)=0).
This problem was discussed here several times, the most recent post is
http://www.mapleprimes.com/posts/207661-Isolation-Of-Sides-Of-The-Surface-On-The-Graph

S is supposed to be the boundary of a convex body having (x0,y0,z0) an interior point and contained in a ball of radius R centered at (x0,y0,z0).
Actually, the procedure also works if the body is only star-shaped with respect to the interior point, and it is also possible to plot only a part of the surface
inside a solid angle centered at (x0,y0,z0).

Usage:
Par3d(f, x=x0, y=y0, z=z0, R, m, n,  theta1 .. theta2,  phi1 .. phi2)

f           is an expression depending on the variables x, y, z
x0, y0, z0  are the coordinates of the interior point
R           is the radius of the ball which contains the surface,
m, n        are the numbers of the grid lines which will be generated
The last two parameters are optional and are used when only a part of S will be parametrized.

The procedure Par3d returns a MESH structure M, which can be plotted with PLOT3D(M).

Par3d :=proc(f,x::`=`,y::`=`,z::`=`,R,m,n,th:=0..2*Pi,ph:=0..Pi)
    local A,i,j, rij,fij,Cth,Sth,Cph,Sph, theta,phi, r;
    A:=Array(1..m+1,1..n+1,1..3,datatype=float[8]);
    for i from 0 to m do for j from 0 to n do
      theta:=op(1,th)+i/m*(op(2,th)-op(1,th));
      phi:=op(1,ph)+j/n*(op(2,ph)-op(1,ph));
      Cth:=evalf(cos(theta)); Sth:=evalf(sin(theta));
      Cph:=evalf(cos(phi));   Sph:=evalf(sin(phi));
      fij:= eval(f, [lhs(x)=rhs(x)+r*Sph*Cth, lhs(y)=rhs(y)+r*Sph*Sth, lhs(z)=rhs(z)+r*Cph]);
      rij:=fsolve(fij,r=0..R);  if [rij]::list(numeric) then rij:=min(rij) fi; 
      if [rij]=[] or not(type(rij,numeric)) then print(['i'=i,'j'=j], fij); rij:=undefined fi; 
      A[i+1,j+1,1]:=evalf(rhs(x)+rij*Sph*Cth);
      A[i+1,j+1,2]:=evalf(rhs(y)+rij*Sph*Sth);
      A[i+1,j+1,3]:=evalf(rhs(z)+rij*Cph);
    od;od:
    MESH(A);
end:

The procedure is not optimized, e.g.
- Cth, etc could be Vectors computed outside the loops
- Some small changes to use evalhf.

###### EXAMPLES ######

f1 := x^2+3*y^2+4*z^2 - x*y - 2*y*z - 10:
plots:-implicitplot3d(f1, x=-5..5, y=-5..5, z=-2..2);

M:=Par3d(f1, x=0,y=0,z=0,5,40,40):
PLOT3D(M);

f2 := x^4+y^4+z^4-1:
M:=Par3d(f2, x=0,y=0,z=0,5,40,40):
PLOT3D(M);

M:=Par3d(f2, x=0,y=0,z=0, 5,40,40, 0..Pi, 0 .. Pi/3): #Plot half of the top only
plots:-display(PLOT3D(M), scaling=constrained);

M:=Par3d(f2,      x=0,y=0,z=0, 5,30,30, 0..Pi, 0 .. Pi):
N:=Par3d(f2+0.01, x=0,y=0,z=0, 5,30,30, 0..Pi, 0 .. Pi):
plots:-display(PLOT3D(M), color=red):
plots:-display(PLOT3D(N), color=green):
plots:-display(%,%%, orientation=[-40,65,10]);

f3 := (x^2+y^2-1)^2+(z+sin(x*y+z))^4-120:
plots:-implicitplot3d(f3, x=-4..4,y=-4..4,z=-5..5, numpoints=10000);

Par3d(f3, x=0,y=0,z=0,5, 30,30):
PLOT3D(%);

Note.
The procedure could be used to plot locally around a point (x0,y0,z0)
One may use the spherical coordinates (theta0,phi0) and then call the procedure taking theta0-a .. theta0+a,  phi0-b, .. phi0+b  for the trailing parameters
The spherical coordonates can be computed using:

ThetaPhi :=proc(x,y,z, X,Y,Z)
    local r:=sqrt((X-x)^2+(Y-y)^2+(Z-z)^2);
    ['theta'=arctan(Y-y,X-x), 'phi'=arccos((Z-z)/r)]
end:

ThetaPhi(10,20,30, 11,21,28);evalf(%);

Download !strange-bug-int.mw