Mathematica, 43 bytes

I've currently got three different solutions at this byte count:

Nest[#+1//.x_/;!(x^2∣(7^x-1)):>x+1&,0,#]&

Nest[#+1//.x_/;Mod[7^x-1,x^2]>0:>x+1&,0,#]&

Nest[#+1//.x_:>x+Sign@Mod[7^x-1,x^2]&,0,#]&

Tcl, 66 bytes

while 1 {if {1>(7**[incr i]-1)%$i**2&&[incr j]==$n} break}
puts $i

Try it online!

PARI/GP, 42 bytes

Pretty straightforward. 1-indexed, though this could easily be changed.

n->=k=1;while(n--,while((7^k++-1)%k^2,));k

or

n->=k=1;for(i=2,n,while((7^k++-1)%k^2,));k

Wonder, 28 bytes

@:^#0(!>@! % - ^7#0 1^#0 2)N

Zero-indexed. Usage:

(@:^#0(!>@! % - ^7#0 1^#0 2)N)2

Filters from the list of natural numbers with a predicate that determines whether x^2 is divisible by 7^x-1, then gets the nth item in that list.

QBIC, 39 bytes

:{~(7^q-1)%(q^2)=0|b=b+1]~b=a|_Xq\q=q+1

I couldn't get it to run in QBasic 4.5, but it seems to run fine in QB64. For some inexplicable reason, QBasic refuses to cleanly divide 13,841,287,200 by 144, but instead gives a remainder of -128. It then returns 16 as the 7th term of this sequence instead of 12...

:{      get N from the command line, start an infinite DO-loop
~       IF
(7^q-1) Part 1 of the formula (Note that 'q' is set to 1 as QBIC starts)
%       Modulus
(q^2)   The second part
=0      has no remainder
|b=b+1  Then, register a 'hit'
]       END IF
~b=a    If we have scored N hits
|_Xq    Quit, printing the last used number (q)
\q=q+1  Else, increase q by 1. 
        [DO-loop and last IF are implicitly closed by QBIC]

Python 2, 57 53 51 bytes

^{-4 bytes thanks to ETHproductions
-2 bytes thanks to TuukkaX}

i=0
g=input()
while g:i+=1;g-=~-7**i%i**2<1
print i

Try it online!
the sequence is 1-indexed

Perl 6, 35 34 bytes

{grep({(7**$_-1)%%$_²},^∞)[$_]}

0-indexed.

Shaved off one byte thanks to Brad Gilbert.

Python 2, 48 46 bytes

Thanks to @Dennis for -2 bytes!

f=lambda n,i=1:n and-~f(n-(~-7**i%i**2<1),i+1)

A one-indexed recursive function that takes input via argument and returns the result.

Try it online! (Recursion limit increased to allow the final test case to run)

How it works

n is the desired index, and i is the counting variable.

The expression ~-7**i%i**2<1 returns True (equivalent to 1) if i^2 divides 7^i - 1, and False (equivalent to 0) otherwise. Each time the function is called, the result of the expression is subtracted from n, decrementing n each time a hit is found; i is also incremented.

The short-circuiting behaviour of and means that when n is 0, 0 is returned; this is the base case. Once this is reached, recursion stops, and the current value of i is returned by the original function call. Rather than explicitly using i, this is done using the fact that for each function call, an increment has been performed using the -~ in front of the call; incrementing 0 i times gives i, as required.

Powershell, too many bytes

Just to see if it was possible and it is.

[System.Linq.Enumerable]::Range(1,10000)|?{[System.Numerics.BigInteger]::Remainder([System.Numerics.BigInteger]::Pow(7,$_)-1,$_*$_) -eq 0}

Perl 5, 35 Bytes

Well, this was missing, so here it is:

map{$_ if!((7**$_-1)%($_**2))}1..<>

Clojure, 83 bytes

(fn[n](nth(filter #(= 0(rem(-(reduce *(repeat % 7N))1)(* % %)))(iterate inc 1N))n))

Try it online!

This builds an infinite list of Java BigIntegers starting at 1 and filters them by the definition. It uses zero-based indexing to select the n^th value from the filtered list.

JavaScript (ES7), 40 bytes

f=(n,i=1)=>n?f(n-!((7**++i-1)%i**2),i):i

This loses precision fairly quickly due to the fact that JS loses precision by 7**19. Here's a nearly arbitrary-precision ES6 version:

f=(n,i=0)=>n?f(n-!(~-(s=++i*i,g=j=>j?g(j-1)*7%s:1)(i)%s),i):i

This finishes within about a second for test case 31.

A few longer approaches:

f=(n,i=0)=>n?f(n-!(~-(s=>g=j=>j?g(j-1)*7%s:1)(++i*i)(i)%s),i):i
f=(n,i=0)=>n?f(n-!(s=++i*i,g=(j,n=1)=>j?g(j-1,n*7%s):~-n%s)(i),i):i
f=(n,i=0)=>n?f(n-!(s=>g=(j,n=1)=>j?g(j-1,n*7%s):~-n%s)(++i*i)(i),i):i

Pyke, 10 bytes

~1IX7i^tR%

Try it here!

~1         -  infinite(1)
  IX7i^tR% - filter(^, not V)
    7i^    -    7**i
       t   -   ^-1
        R% -  ^ % v
   X       -   i**2
           - implicit: print nth value

R, 35 bytes

This only works for n<=8.

z=1:20;which(!(7^z-1)%%z^2)[scan()]

However, here's a longer version which works for n<=25, for 50 bytes:

z=1:1e6;which(gmp::as.bigz(7^z-1)%%z^2==0)[scan()]

PHP, 47 49 bytes

while($n<$argv[1])$n+=(7**++$x-1)%$x**2<1;echo$x;

Only works for n<9 (7**9 is larger than PHP_INT_MAX with 64bit)

62 bytes using arbitrary length integers: (not tested; PHP on my machine doesn´t have bcmath)

for($x=$n=1;$n<$argv[1];)$n+=bcpowmod(7,++$x,$x**2)==1;echo$x;

Run with php -nr '<code>' <n>.

pseudo code

implicit: $x = 0, $n = 0
while $n < first command line argument
    increment $x
    if equation is satisfied
        increment $n
print $x

Python 3, 45 bytes

f=lambda n,k=2:n<2or-~f(n-(7**k%k**2==1),k+1)

Return True for input 1, which is allowed by default.

Try it online!

Python 2, 57 bytes

This takes a really, really long time for large values. It also uses plenty of memory, because it builds the entire list way farther than necessary. The result is zero-indexed.

lambda n:[x for x in range(1,2**n+1)if(7**x-1)%x**2<1][n]

Try it online

05AB1E, 11 bytes

µ7Nm<NnÖiN¼

Try it online!

For some reason I can't get ½ to work in µ7Nm<NnÖ½N or I'd be tied with Pyth.

µ           # Loop until the counter equals n.
 7Nm<       # Calc 7^x+1.
     Nn     # Calc x^2.
       Ö    # Check divisibility.
        iN¼ # If divisible, push current x and increment counter.
            # Implicit loop end.
            # Implicitly return top of stack (x)

.

Pyth, 10 bytes

e.f!%t^7Z*

A program that takes input of an integer and prints a one-indexed value.

Try it online!

How it works

e.f!%t^7Z*     Program. Input: Q
e.f!%t^7Z*ZZQ  Implicit variable fill
               Implicitly print
e              the last
 .f         Q  of the first Q positive integers Z
     t^7Z      for which 7^Z - 1
    %          mod
         *ZZ   Z^2
   !           is zero

Haskell, 34 bytes

([x|x<-[1..],mod(7^x-1)(x^2)<1]!!)

This uses 0-based indexing. Usage example: ([x|x<-[1..],mod(7^x-1)(x^2)<1]!!) 30 -> 1140.

It's a direct implementation of the definition. It builds a list of all numbers x and picks the nth.