You can use the codegen package, g1:=codegen[makeproc](ff, [x]) and g2:=prep2trans(g1)

Then plot(g2, x0-3 .. x0+3) works for me

treat each term and change variables

Download MP_236350-How-To-Evaluate-A-Onedimensional-Improper.mw

For what its worth and even if this is not an Excel forum: here a sheet, zipped.

The sheet contains a function "KeplerZero" which solve the equation F(x,c) = x - c*sin(x) = 0 for a given c, 1 < c, 0 < x < Pi.

At the boundaries one can use a Taylor series for F, it is alternating and therefore that approximation error can be estimated.

In x = 0 one takes series = polynom of degree 5, and the task easily reduced to solve a quadratic equation.

In x = Pi one takes degree 3, a depressed cubic with 3 real roots (Descartes rule of sign).

Between those ranges one uses a root finder, already a bisection method will do

Note: for increasing c the task runs into catastrophic cancellations since both summands are necessarily of same size
and achieving zero means that their leading digits have to cancel out. Meaning that the true value can not be seen.

PS: one can use it with Open Office / Libre Office as well, the code should work.

axelvogt.de/temp/KeplerZero_BiSection.xls (I can not upload it here ...)

circle([0,1], 2, color=blue), I presume

If it is true then you can just define it

It works. Note from the help: "If X is a list, then Shuffle returns a new copy either way."

Hence you want to use the shuffeled list, not the original one.

Download MP_236210.mw

For example using fsolve( cos(s)*cosh(s) - 1, s= 0 .. 10, maxsols=infinity );

May be that you have to refine those solutions.

For large values you may re-scale, say deviding by 1/2*exp(s), which has the same zeros - but does not "explode" as it simplifies to cos(s)+cos(s)*exp(-2*s)-2*exp(-s)

map( lhs, eq1 )

Pn is already defined, as a rational function, having a possible pole of order 1 in +- I.

Computing the residue in that points it turns out to be zero. So it is a polynomial.

If you write simplify(P) then you will see that it has a factor (x+y+1)^3. Done

I changed variables and equations by x6,x7,x8:=x,y,z and then by y=x+r, z=y+t=z=x+r+t and conditions 0<x, 0<r, 0<t, x+r+t<1.

Plotting indicates that these surfaces to not intersect in the desired range. This is not a prove and those plots will not catch isolated solutions.

MP_235991_Intersect_equations.mw

For example

singular(1/sin(x)/cos(x),x) or singular(1/sin(x+cos(Pi*x)),x);

n:='n';
integrand:=(2*x^n+1)/(x^(n+1)+x);
int(integrand,x); simplify(%) assuming 0<n;

giving ln(x) + ln(x^n + 1)/n

Do you really mean that or numerical zero (may be after "fnormal") ?

I get a solution -ln(inf^3*r + (6*I*r*t - r^3 - 2*eps^3)*inf + 2*eps^3*r) +
ln(-(-inf + r)*(-2*eps^3 + inf^2*r + inf*r^2)*exp(exp(-r^2/2)/4) + 6*I*inf*r*t)

which would do it numerical feeding your values (after "fnormal"), but that
boils down to exp(-r^2/2)/4 for eps --> 0, inf --> infinity

You can try gfun, an external package

with(gfun):
-4*RootOf(_Z^3 - 3*_Z^2 - 10*_Z - 1)^2 +
19*RootOf(_Z^3 - 3*_Z^2 - 10*_Z - 1) + 3:
algfuntoalgeq(%, y(x));

                     3       2             
                    y  + 50 y  + 75 y - 125

https://perso.ens-lyon.fr/bruno.salvy/software/the-gfun-package/