Let  a  be the number of points on one dice, and   b  be the number of points on another one.
 

Download LE.mw

You can use an inert form for subtraction to keep the expression unchanged. For any calculations and transformations use the  value  command:

restart;
expr:=arccos(p%-a);
value(expr);
value(eval(expr,[p=0,a=1]));

_B1  is a binary variable that is 0 or 1 :

restart;
ode:=diff(y(x),x)-y(x)/x+csc(y(x)/x)=0;
sol:=dsolve([ode,y(1)=0]);
about(_B1);
sol1:=eval(sol,_B1=0);
sol2:=eval(sol,_B1=1);
odetest(sol1,ode);
odetest(sol2,ode);

Here is a solution in Maple:
 

Download nlsystem.mw

restart;
plots:-animate(plot,[cos(2*theta),theta=0..a, coords=polar], a=0..2*Pi);

We can significantly increase efficiency if we use symmetries: from one solution we can get  6 * 6 = 36 solutions by permutations of the first three and last three numbers. Instead of loops, it’s slightly more efficient and more compact to use nested seq . So we get the unique solution:

Download isys.mw

Your surface is a torus:

restart;
plot3d(eval([x,sqrt(x^2+y^2)*cos(s),sqrt(x^2+y^2)*sin(s)],[x=cos(t),y=2+sin(t)]), t=0..2*Pi, s=0..2*Pi, scaling=constrained, labels=[x,y,z]);

Download shade.mw

Here is another shorter way:
 

Download shade1.mw

Here is a more automatic plotting for the intersection:

restart;
r1:=2:  r2:=2*(1-cos(theta)):
plot(min(r1,r2), theta=0..2*Pi, color=green, coords=polar, filled, scaling=constrained);

Edit.

The following approach is sometimes useful:

restart;
eq := x^2+floor(x)-10:
Student:-Calculus1:-Roots(eq);
identify(%);

The animation you see here is saved in GIF format. Each polygon from the list remains on the display for exactly 1 second (10 frames per second by default):

restart;
OneFrame:=proc(n)
local Circle, Polygon, Text, S;
uses plottools, plots;
Circle:=plot([cos(t),sin(t),t=0..2*Pi], color=grey, thickness=2);
Polygon:=polygon([seq([cos(2*Pi*k/n),sin(2*Pi*k/n)], k=1..n)], color=yellow);
S:=evalf(n*sin(2*Pi/n)/2);
Text:=textplot([[-0.7,-1.1,"Area of Circle = 3.14"],[0,-1.1,typeset("n","=",n)],[0.7,-1.1,typeset("Area of Polygon","=",evalf[3](S))]], font=[times,16]);

display(Text,Circle,Polygon, scaling=constrained, axes=none, title="Archimedes Approximation for Pi \n", titlefont=[times,bold,18]);
end proc:

plots:-display(seq(OneFrame(n)$10, n=[3,6,9,12,15,18,21,24,27,30,60]), insequence, size=[600,600]); 

Download diffeq.mw

Maple has the command  Student:-LinearAlgebra:-LinearSolveTutor  that solves systems of linear equations step by step, but unfortunately only if the matrix of the system is no more than 5 by 5. Below is a step-by-step solution using the Jordan Gauss method using my program  JordanGausse (all comments are in Russian):

system.mw

reduce the range along the vertical axis (y-range) while maintaining  x-range. Also use a list rather than a set when specifying functions so that all options match the corresponding function:

restart;
plot([0, 2*x^2, 2*x^2 - 2*x^3 + 8/3*x^4 - 4*x^5], x = -10 .. 10, y=-1000..1000, color = ["DarkGreen", "CornflowerBlue", "Burgundy"], axes=box);

Here is another way by the using the  applyrule  command:

restart;
expr:=1/exp(z)*arcsinh(x*exp(C[1]))+x*sin(exp(x))+3*exp(C[1]*y)*sqrt(sin(exp(3*h)));
applyrule(exp(t::anything)=Z, expr);

I think that there are infinitely many ellipses inscribed in a given triangle. But among them there will be only one that touches the sides of the triangle in their midpoints. It is called the Steiner ellipse. See  https://en.wikipedia.org/wiki/Steiner_inellipse

Here is an example of plotting this ellipse:
 

As for orthoptic circle for an ellipse see  https://en.wikipedia.org/wiki/Orthoptic_(geometry)#:~:text=Generalizations%3A,fixed%20angle%20(see%20below).&text=Thales'%20theorem%20on%20a%20chord,two%20points%20P%20and%20Q.

Download Steiner_ellipse.mw