@sursumCorda 
I used Maple 2018 on an old Windows tablet  and I obtain 
limit(1/f, x=infinity) = 0.

@petit loup   "#" starts a comment, so actually method='integer'  is not needed when the matrix has integer entries.
For our matrix A,
S,U,V := SmithForm(A, output=['S,U,V'], method='integer');
and
S,U,V := SmithForm(A, output=['S,U,V']);
give the same result.
 

@mmcdara Actually, when limit returns an interval a..b, this means that the limit does not exist; mathematically this is equivalent to the undefined result, but Maple also says that the limit points (aka accumulation points) are contained in the interval a..b.

Note that e.g. 

limit(3*sin(1/x)-4*cos(1/x), x=0);
    -7 .. 7

but actually the set of accumulation points is the interval -5 .. 5, which Maple cannot find.

It would be nice to have in Maple the comands  sup, inf, limsup, liminf for functions and sequences, just like in Mathematica!

@mmcdara Probably the designers consider that <...> is enough, e.g. the first method or

<M[1,1],M[1,2]; M[2,1],M[2,2]>;

@Mikey_W 

restart;
m:=1000: n:=1000: A:=LinearAlgebra:-RandomMatrix(m, n):
x:=Pi/77:
    
f:=proc(x,A) local i,j,s:=0; 
for i to m do for j to n do
  s:=s+cos(sin(x+A[i,j]));
od od;
s
end proc:

CodeTools:-Usage( evalf(f(x,A)) );
## memory used=6.11GiB, alloc change=33.00MiB, cpu time=10.89s, real time=10.90s, gc time=1.70s

##                        7.639612369 10^5 

CodeTools:-Usage( evalhf(f(x,A)) );
## memory used=7.63MiB, alloc change=7.63MiB, cpu time=265.00ms, real time=266.00ms, gc time=0ns

##                    7.63961236894650501 10^5 

@acer  No, please do.

@sursumCorda Corrected and .mw inserted.

Download sysok_vv.mw

@mmcdara If you don't know an IC for some x>1, you cannot find a solution.
For example, 
sol3 := piecewise( x<1, 1/(1-x), 1/(a-x) )

is a solution for any a<1.

Removing the conditions psi1(1)=0, psi2(1)=0, the solution is unique in terms of e1, e2 and can be obtaind symbolically.

Now, imposing psi1(1)=0, psi2(1)=0 it results a system in {e1, e2} which can be solved numerically.

@Carl Love I do not have a solid theoretical justification. The vector <x,y> can ve represented (via isolve) as p*v1 + q*v2; p,q in Z, v1, v2 in Z^2.  For our problem, it will be convenient to have small p,q when x and y are small. This happens when v1, v2 are replaced with the simplified equivalent basis given by LLL.

I think that indeed R_max could be estimated using the properties of the LLL basis, but I did not (and probably will not) try to investigate this.

@Carl Love Nice idea, I have used it in my first version but without chrem.
Unfortunately, for larger problems it is very slow.
For example, I could not wait MinY for
eqs := [154*x+69*y, 13*x+716*y, 23*x+3059*y, 2295*x + 4522*y, 6479*x + 5396*y, x+2*y]:
M:= [7^3, 13^3, 23^3, 17^3, 29^3, 31^3]

but Miny finished in less the 1 second (for R=100) or 5 seconds for R=1000 (default).
 

Congratulations for  mapletransactions and for your interesting contribution on Special Functions.

P.S. Concerning Tom Rocks Maths, I wonder whether such "performances" are really necessary to attract young math students.

@Jesús Guillera  I guess your comment is instead of a "thanks".
Have you tried my proc? It works correctly.
In my second example, the last modulus is 29^3. (For your 19^3 ==> x = 1103, y = 26390).
 

@Jesús Guillera  Without a minimization the solution is not unique (modulo the product of the moduli) and you already have here an example.