linear Diophantine equations

When re-expressed in maths terms, the current Riddler is about finding a sequence x⁰,x¹,…,x⁷ of integers such that

x⁰=7x¹+1
6x¹=7x²+1
…
6x⁶=7x⁷+1
6x⁷=7x⁸

which turns into a linear equation with integer valued solutions, or a system of linear Diophantine equation. Which can be easily solved by brute-force R coding:

A=matrix(0,7,7)
for (i in 1:7) A[i,i]=6
for (i in 1:6) A[i,i+1]=-7
for (x in 1:1e6){
  zol=solve(a=A,b=c(rep(1,6),7*x))
  if (max(abs(zol-round(zol)))<1e-3) print(x)}
x=39990 #x8=5.6.31.43
7*solve(a=A,b=c(rep(1,6),7*x))[1]+1 #x0

which produces x⁰=823537. But it would be nicer to directly solve the linear system under the constraint. For instance, the inverse of the matrix A above is an upper triangular matrix with (upper-)diagonals

1/6, 7/6², 7²/6³,…,7⁶/6⁷

but this does not help considerably, except for x⁸ to be solutions to 7 equations involving powers of 6 and 7… This system of equations can be solved by successive substitutions but this still feels very pedestrian!