@Carl Love Oops, I just found it and edited my reply. Thanks.

@Carl Love I believe that this is related. 

restart;
libname := subs(FileTools:-JoinPath([kernelopts(':-toolboxdir' = "Physics Updates"), "lib"]) = NULL, [(tmp := libname)])[]:
inttrans[':-invlaplace'](1/(s + a), s, t);
                           exp(-a t)

libname := tmp:
forget(inttrans[':-invlaplace']);
inttrans[':-invlaplace'](1/(s + a), s, t);
                            exp(a t)

In addition, 

> inttrans['invlaplace'](inttrans['laplace'](ln(t^2+1), t, s), s, t); # Maple 2023.1 
 = 
                              / 2    \
                          2 ln\t  + 1/

So Maple thinks that ㏑(t²＋1)＝2㏑(t² + 1)…? 

Additional examples: 
 

Download ilaplace.mw

(1) is correct, but (2), (3), (4), and (5) are wrong.
And for uknown reason, Maple cannot calculate (6), (7), (9), and (10).

@Thomas Richard Right. But consequently, if one has installed Python before, one will have to add a multiple Python on the machine (when installing Maple). Will such a Python installation be optional in future Maple releases?

I think that you mean: 

Python:-ImportModule("import math, numpy as np");
Python:-DefineFunction("def Concentration_calculation(C_0, Q, V_r, m_b, rho, R, Gamma_i, delta_t=1):\n\tt = np.arange(0, 360*60, delta_t)\n\tC_i = [C_0]\n\tfor i in range(len(t)-1):\n\t\tdC = -(Q/V_r)*(1-math.exp(-3*m_b*math.sqrt(Gamma_i/(rho*t[i+1]))/(math.sqrt(math.pi)*Q*R)))*C_i[i]\n\t\tC_i.append(C_i[i] + dC*delta_t)\n\treturn t, C_i"):

@Carl Love Many thanks. This is an impressive procedure. (At least the computation can be done in one-third of a minute!) 

The fourth parameter `D` being 0 appears to mean that not all digits can be used, but there exists a bug in your program:

A161786__2(1, 1, 10, 0); # though it is not hard to fix it 
Error, (in A161786__2) Array index out of range

.

@dharr Thanks. In practice, I find that most use of “convert(p, 'base', 10)” can be simply replaced by faster “[`convert/base`:-MakeSplit(length(p), 1, 10)(p)]”, but to use it, one has to set "kernelopts('opaquemodules' = false):" in advance (which is somewhat dangerous). Besides, after converting integers to strings, I also try to call StringTools:-RegMatch directly, but Maple's regular expression engine seems to be not optimized enough ("Regular expressions" are extremely arcane, but they are known for being very powerful as well.). 
As for the `ithprime`, I believe that its inefficiency is because it's "implemented in the high-level Maple programming language". But I don't know why such a fundamental function (for a math software) is not built into the core of the Maple computation engine …. 

@mmcdara Thanks. I think that your `B1` and `B2` can be obtained more directly (using NumberTheory:-PrimeCounting(N1) + 1 and NumberTheory:-PrimeCounting(10*N2 - 1)).
But in this question, the desired interval is p_1000000～p_1999999 (corresponding to B1 = 10^6 and B2 = 2*10^6 - 1). (I consider testing this range because it may be suited for a benchmark problem—neither too large nor too small.) Anyway, though I ask for a large-scale problem, I'm happy to find that your code performs better on small numbers.

How about 

RealDomain:-solve({x + abs(y + 1) = 1, y + abs(z + 2) = 1, z + abs(x - 2) = 1});
       {x = -abs(y+1)+1, z = -abs(y+1), -1 <= y, y <= 1}

@vv For example, 

But one can easily check that (x, y)＝(78917749904961, 78919591651475) is better. Regrettably, I have no time to wait for a larger R. (Note that minimizing x＋y will lead to a bigger solution (x, y) = (158339291, 112245157230352).)

@Carl Love Here is the Maple code: 

mfoo := proc(` $`)::posint:
	local bound::posint, s::object := Session('logic' = "QF_LIA");
	uses SMTLIB;
	s:-Assert(And((154*x + 69*y) - 0 = 7**3*m__0, (13*x + 716*y) - 0 = 13**3*m__1, (23*x + 3059*y) - 0 = 23**3*m__2, x > 0, x < y));
	for bound do
		s:-Push();
		s:-Assert(y = bound);
		ifelse(s:-Satisfiable(), RETURN(bound), s:-Pop())
	od
end:

CodeTools:-Usage(mfoo(), iterations = 1);
 = 
memory used=0.97GiB, alloc change=0 bytes, cpu time=88.03s, real time=88.82s, gc time=2.95s

                             18837

And a practically word-for-word translation: 

Python:-Verbose((*not *)false);
ssystem(sprintf("PowerShell -Command \"& '%spython.exe' -m pip install -U --force-reinstall z3-solver --no-cache-dir\"", StringTools:-Subs("bin" = "Python", kernelopts(bindir)))): # Remember to uninstall `z3-solver` later! 
Python:-Import("from z3 import *"):
Python:-DefineFunction("def pfoo() -> Int:\n\tx, y = Ints('x y'); m = IntVector('m', 3)\n\ts = SolverFor(\"QF_LIA\")\n\ts.insert((154*x + 69*y) - 0 == 7**3*m[0], (13*x + 716*y) - 0 == 13**3*m[1], (23*x + 3059*y) - 0 == 23**3*m[2], x > 0, x < y)\n\tbound = 0\n\twhile True:\n\t\ts.push()\n\t\ts.insert(y == (bound := bound + 1))\n\t\tif s.check() != sat:\n\t\t\ts.pop()\n\t\telse:\n\t\t\treturn bound"):
CodeTools:-Usage(Python:-EvalFunction("pfoo"), 'iterations' = 1);
memory used=1.20KiB, alloc change=0 bytes, cpu time=31.00ms, real time=8.60s, gc time=0ns

                             18837

There exist some other state-of-the-art SMT solvers (e.g., Yices 2, Boolector, Bitwuzla, and cvc5); nevertheless, the first three are not convenient to use on Windows (whereas the last one provides pre-built binaries for Windows). If I have time, I shall test cvc5 (whose input can be obtained by SMTLIB:-ToString) in Maple. 

@Carl Love Thanks for your advice.
As I see it, the program should not be too "special" since if the OP plans to additionally minimize x＋y instead of the current y (Note that their results are different. One is (x, y) = (15893, 18837); the other is (x, y) = (1103, 26390).), the objective expression will no longer be a multiple of 7×13 (though some other tricks will remain available). So I attempt to concentrate on defining OP's central task and write code that captures concepts rather than burning in a particular generation of specific algorithms. 

Anyway, I believe that there is something wrong with Maple's SMTLIB package. According to 

showstat(SMTLIB:-Execute, 1);

SMTLIB:-Execute := proc(smt::string, {getassignment::truefalse := false, getmodel::truefalse := false, getproof::truefalse := false, getunsatcore::truefalse := false, getvalue::sequential(name) := [], solver::Or(identical(maple,z3),string) := ':-z3', z3options::sequential(equation) := []})
   1   if solver = (':-z3') then
           ...
       elif solver = (':-maple') then
           ...
       else
           ...
       end if
end proc

Maple by default makes use of the external Z3 prover. However, even the naïve Python analogue runs considerably faster than the raw Maple implementation: 

What happened in Maple?! Maybe certain overheads of linking and calling external routines become critical here. (Recall that Python's native loops are known for their poor performance!)

@Axel Vogt I'm not sure what happened. On my computer: 
 

Download question237110.mw

This is far slower than @vv's one. (Why?)

@vv If I have understood right, (x, y) = (15893, 18837) should be more "minimal" than (x, y) = (1103, 26390). (Note that OP's definition is: x and y positive, x＜y, and y as small as possible.) However, computing with finite precision might lead to wrong results.

For medium numbers of nodes, one may use the simulated annealing or the genetic algorithm to find a route whose length is at least close to the minimum as well (see Approximate a Solution to Traveling Salesman Problem). But strangely, though a number of built-in procedures use compiled subroutines from NAG Library as their underlying engine and to carry out and accelerate computations and there does exist a function h03bbc provided by the NAG that approximately solves the symmetric TSP, GraphTheory:-TravelingSalesman still does not make use of this NAG solver. (Why? I read that "the incorporation of the NAG library into Maple is an ongoing project" (cf. work on including NAG routines for other types of numeric computation is continuing); maybe the developers of GraphTheory package are not very familiar with those of the LinearAlgebra package ….) 

Edit. Here is an additional comment. MatLab's marketing literature claims that classical TSPs can also be solved by mapping them to the so-called Quadratic Unconstrained Binary Optimization (QUBO, not MILP) problems, which appear to be a terra incognita of Maple's Optimization package.