I wonder, if there are examples of Diophantine equations (or systems of such equations) with integer coefficients fitting on a few lines that have been proven to have a finite, but really huge number of solutions?
Are there ones with so large number of solutions that we cannot write any explicit upper bound for this number using Conway chained arrow notation?
Update: I am also interested in equations with few solutions but where a value in a solution is very large itself.
On
There are many such equations. Here is one such equation.
Note that the number of ways of writing a number $n$ as a sum of $2$ squares, where $$n=2^{a_0}\left(p_1^{a_1}p_2^{a_2}\ldots p_k^{a_k} \right)^2 q_1^{b_1} \cdots q_s^{b_s}$$ $p_i \equiv 3\pmod4$, $q_j \equiv 1 \pmod4$ and $p_i,q_j$ are primes, is given by $$f(n) = 4B$$ where $B = (b_1+1)(b_2+1)\cdots(b_s+1)$. Choose a large enough $B$. Construct a number $n$ with the given $B$ (You have plenty of options of creating $n$). Form the diophantine equation $x^2+y^2 = n$. This has $4B$ solutions.
You can create similar examples like $x^2 + 2y^2 = n$, $x^2 + 3y^2 = n$ and letting $n$ to be large with some constraints, you can get as many number of integer solutions, you want.
On
I think the rough answer to your question is that "there are no such natural equations." Let me try to justify this heuristic.
Caporaso, Harris, Mazur; and Baker's method: Consider diophantine equations corresponding to curves, for example equations of the form:
$$C: f(x,y) = 0$$
where $f(x,y)$ is a polynomial with integer coefficients. One knows that the number of solutions (with $x$ and $y$ in $\mathbf{Z}$ or $\mathbf{Q}$) is --- to some extent --- determined by the complex geometry of $C$. If the genus $g$ of $C$ is greater than two, then Faltings proved that $C$ has only finitely many rational points. However, much more (may) be true. Assuming a conjecture of Lang, it was proved by Caporaso, Harris, and Mazur (J. Amer. Math. Soc. 10 (1997), 1-35) that the number of rational points $\#C(\mathbf{Q})$ is bounded by a function $A(g)$ which only depends on $g$, not on $C$. This doesn't prove that $A(2)$ [for example] is not enormous, but all known lower bounds in $g \ge 2$ on the number of rational points are at most polynomial. For example, I think the largest number of rational points on a genus two curve that has ever been found has at most a few hundred rational points. Since the genus is controlled by the degree of $f$, this suggests that it is highly unlikely to find an equation with small coefficients and degree which has absolutely massive solutions, even if one makes the coefficients big. Note that this is only a heuristic, but it strongly suggests that there are no such equations. It certainly says that no such equations are currently known.
One is in slightly better shape (or worse, depending on what you're trying to do) if one restricts to integral points of $C$. For certain classes of equations (say $F(x,y) = m$ for homogenous $F$) there are explicit bounds on the number of integral solutions in terms of the coefficients coming from Baker's theorem (on linear forms in logarithms). The bounds are not great (perhaps super-super exponential in the coefficients), but they certainly preclude numbers of the size you are interested arising from anything that one can sensibly write down. Moreover, one conjecturally expects that Baker's bounds are not optimal.
For genus 1 and 0, the situation is similar. There's difference here again between whether one wants to consider only integral solutions or rational solutions, but the result (in either case) is that if one insists that there are only finitely many solutions (in either integers or rationals), then there is a bound for the largest such solution in terms of the coefficients (which, using Baker's method for integral points, might be quite large, but is still tiny compared to the numbers you are discussing).
Moral: If one assumes Lang's conjecture (and we certainly don't have any evidence that it is false), there is still heuristic evidence to suggest that the the types of equations you are looking for do not exist, either for curves or higher dimensional varieties: this is very speculative, but certainly means that writing down any such equations is either impossible or very hard.
Exponential Diophantine Equations: There are other flavours of Diophantine equations besides algebraic varieties. For example, you could consider exponential Diophantine equations, e.g.: $2^n = x^2 + 23$. In this case, one can often convert solutions to these equations into subsets of polynomial equations of genus at least two. For example, the previous equation gives solutions to one of the five genus three curves: $$Ay^5 = x^2 + 23, \qquad A = 1,2,4,8,16.$$
Thus one can often "reduce" to the case of curves. Alternatively, Baker's method can often be used directly on exponential equations to give upper bounds.
The phenomena of a big "smallest" solution: There are some equations which have infinitely many solutions of which the "smallest" is quite large, for example, the first non-trivial solution to Pell's equation $x^2 - 61 y^2 = 1$ has $$[x,y] = [1766319049, 226153980].$$ However, in these cases (and in similar cases coming from elliptic curves) one at least expects (and knows for Pell) that there exists a bound on the regulator which is polynomial in the coefficients, where the regulator involves the logarithm of the entries. So this gives upper bounds for primitive solutions which are exponential in a polynomial in the coefficients, still enough to prevent super-huge numbers arising as the smallest interesting solution.
Summary: There are many other flavours of Diophantine equation one can write down, but I think the summary is that the conjectures of Lang suggest, on some heuristic level, that the type of phenomena you are looking for simply does not occur naturally. There may be a way of embedding Ackerman's function in to some system of equations (but not into polynomial equations), but I think you would consider that cheating.
On
The smallest solution to the system $$x^2+101y^2=z^2,\quad x^2-101y^2=t^2$$ is $$x=2,015,242,462,949,760,001,961;\quad y=118,171,431,852,779,451,900;$$ $$z=2,339,148,435,306,225,006,961;\quad t=1,628,124,370,727,269,996,961$$ according to Wolfram.
On
While I largely agree with @Qbert's sentiment that you won't find anything that's hard to bound with Conway chained arrow notation (note that Ackermann's function and Graham's number only need 2 or 3 arrows, so even allowing the shortcut "where G is Graham's number" won't get you to the next line), here is something explicit that at least you might want to use two arrows for. $$ a^3-k^2=-11492\\ x^2-(a^2-1)y^2=1\\ u^2-(a^2-1)v^2=1\\ s^2-(b^2-1)t^2=1\\ v=ry^2 \quad b=1+4py \quad b=a+qu \\ s=x+cu \quad y=k+e-1 \\ t=k+4(d-1)y\\ x = w+z \\ f^3-g^2=w $$
I consider solutions in positive integers. (You can add constraints like $ w=1+A^2+B^2+C^2+D^2 $ where necessary to force variables to be positive if you prefer to work over all of $\mathbb Z$.)
The first equation has a solution with $a=154319269$ as given by Elkies here (and also others, e.g. $a=13,k=117$). The next six lines are exactly the system given by Martin Davis for Theorem 3.1 in "Hilbert's Tenth Problem is Unsolved" and only have one solution for given $a,k$ when $x$ is exactly the $k$th smallest solution to $x^2-(a^2-1)y^2=1$, in which case $x\ge a^k$. The next lets $w$ range over all $1\le w < x$.
The first and last equations are Mordell curves and have only finitely many solutions for $w$ fixed. Taking the value for $a$ given above, there is a solution with $x>10^{10^{13}}$. Then I couldn't really say how many solutions there might be to the last equation, but the largest is at least $\lfloor\sqrt[3]{x-2}\rfloor$, and based on Marshall Hall's conjecture we might expect it to be greater than $x^6$.
But in any case a bound from Baker given in section V.8 of Silverman & Tate would require $$ g\le\exp((10^6w)^{10^6})<3^{10^{10^{20}}}<10\rightarrow 5 \rightarrow 2 $$ using Conway's notation (and assuming the $a$ given above is the largest solution to the first equation).
I couldn't give you an explicit equation, but if you're willing to bend the rules a bit, I'm fairly certain that one should exist. Apply Matiyasevich's theorem to the busy beaver problem: for any $n$, the set of running times of all halting binary Turing machines with $n$ states is finite and Diophantine, and therefore represented by the positive values taken by some integer polynomial $P(x_1,\ldots,x_k)$.
Then one could choose the Diophantine equation $P(x_1,\ldots,x_k) = 1+a^2+b^2+c^2+d^2$, with the proviso that there could be infinitely many solutions (though perhaps someone familiar with the MRDP construction could correct me), but there are only finitely many distinct values of $(a,b,c,d)$. And the busy beaver function grows absurdly (non-computably) fast, so it shouldn't take very long for the values of $a$, $b$, $c$, $d$ to outstrip anything that arrow notation can describe.