APL(NARS), 48 chars

r←h w;a
a←⍬⋄r←1
→2×⍳∨/0=a∣r+←4
a,←r⋄→2×⍳¯1<w-←1

// +/8 8 15 17=48

  )box on
Was OFF
  0=⍬∣5
┌0─┐
│ 0│
└~─┘
  ∨/0=⍬∣5
0
~

this means that exit from first loop in line 2, and enter in line 3 in the second loop with r=5

  h¨0..10
┌11─────────────────────────────┐
│ 5 9 13 17 21 29 33 37 41 49 53│
└~──────────────────────────────┘
  h 500
3749
~~~~
  h 10000
96589
~~~~~

Vyxal 3, 13 bytes

5⎄{:K#¤∩|4+]i

Vyxal It Online!

5⎄{:K#¤∩|4+]i­⁡​‎‎⁡⁠⁤⁡‏‏​⁡⁠⁡‌⁢​‎‎⁡⁠⁢‏⁠‎⁡⁠⁣⁤‏‏​⁡⁠⁡‌⁣​‎‎⁡⁠⁡‏‏​⁡⁠⁡‌⁤​‎‎⁡⁠⁣‏⁠⁠‎⁡⁠⁣⁤‏‏​⁡⁠⁡‌⁢⁡​‎‎⁡⁠⁢⁤‏‏​⁡⁠⁡‌⁢⁢​‎‎⁡⁠⁢⁢‏⁠‎⁡⁠⁢⁣‏‏​⁡⁠⁡‌⁢⁣​‎‎⁡⁠⁤‏⁠‎⁡⁠⁢⁡‏‏​⁡⁠⁡‌⁢⁤​‎‎⁡⁠⁣⁡‏‏​⁡⁠⁡‌⁣⁡​‎‎⁡⁠⁣⁢‏⁠‎⁡⁠⁣⁣‏‏​⁡⁠⁡‌­
             i # ‎⁡index implicit input into
 ⎄         ]  # ‎⁢generator
5              # ‎⁣with initial vector 5
  {        ]   # ‎⁤while:
       ∩       # ‎⁢⁡the intersection of
     #¤        # ‎⁢⁢all previosuly generated
   :K          # ‎⁢⁣and the number's factors
        |      # ‎⁢⁤is non-empty
         4+    # ‎⁣⁡add 4 to the number
💎

Created with the help of Luminespire.

<script type="vyxal3">
5⎄{:K#¤∩|4+]i
</script>
<script>
    args=[["0"],["1"],["2"],["60"]]
</script>
<script src="https://themoonisacheese.github.io/snippeterpreter/snippet.js" type="module"/>

TI-BASIC (TI-84 Plus CE Python), 61 46 bytes

Input G
For(I,1,ᴇ5
If 2=sum(not(remainder(4I+1,seq(4N+1,N,0,I
DS<(G,1
End
1+4I

gets crazily slower when you up the number

AWK, 93 bytes

{for(i=5;i<$1^2;i+=4)b[++x]=i;for(;++j&&y<$1;b[j]==b[k]&&y++)for(k=0;b[j]%b[++k];);}$0=b[j-1]

Attempt This Online!

{for(i=5;i<$1^2;i+=4)    # input squared for sufficiently long list
b[++x]=i;                # create array of hits
for(;++j                 # traverse array
&&y<$1;                  # until we reach our prime
b[j]==b[k]&&y++)         # if prime, increment y
for(k=0;b[j]%b[++k];);}  # primality test
$0=b[j-1]                # set output

JavaScript, 59 58 54 bytes

0-indexed and based on the (possibly erroneous?) observation that, just as all composite numbers are divisible by at least one prime, all Hilbert composites are also divisible by at least one Hilbert prime.

Here's a limited proof of that using the first 20000 Hilbert numbers.

We handle the edge case of 1 being neither prime nor composite below by hardcoding 5 as our stating point.

a=[x=5];f=n=>a[a.every(y=>x%y,x+=4)?a.push(x):n]||f(n)

Try it online!

Explanation

We initialise an array a=[5], which will hold the Hilbert primes, and an integer x=5, which we'll use as a temporary variable for the Hilbert numbers, outside the main function. Usually this wouldn't be possible because our rules require that everything be reusable by successive function calls, but here on each successive call of the function a[n] will either exist and be immediately returned or we'll resume building the array of primes from where we left off on the last call and continue until it does exist.

f, then, is our recursive function taking input via parameter n. On each iteration we increment x by 4 and then check that, for each existing prime y in a, x%y>0.

If it is then we push the new value of x to a, which returns the updated length of a. We use that to index into a but, as JavaScript uses 0-based indexing, a[a.length] will always return undefined, which is falsey, so we fall through the logical OR to execute a recursive call to f(n) and continue to do so until we encounter a Hilbert composite (i.e., for at least one y in a, x%y===0).

When we do encounter a composite we instead index n into a and, if it exists, return that value, or continue recursing until the next composite if it doesn't.

Raku, 35 bytes

grep {$_%%none 5,9...^$_},(5,9...*)

Try it online!

This is an expression for the lazy, infinite list of Hilbert primes.

• 5, 9 ... * is the infinite list of Hilbert numbers. (Raku infers that it's an arithmetic sequence from the first two elements.)
• grep { $_ %% none 5, 9 ...^ $_ } filters that list to those which are divisible (%%) by none of the Hilbert numbers up to, but not including (...^), the number being tested.

Python, 80 bytes

f=lambda I,i=0,H=[2]:(i>I)*H[0]or f(I,i+all((p:=4*len(H)+1)%h for h in H),[p]+H)

Attempt This Online!

0-indexed (but it would be trivial to change to be 1-indexed). Stores the previous Hilbert numbers in H, and stores how many Hilbert primes have been found so far in i.

Perl 5 -p, 56 bytes

$k=$.+=4;1while($k-=4)>1&&$.%$k;($k>1||--$_)&&redo;$_=$.

Try it online!

Jelly, 17 15 bytes

ŻḤḤ‘ḍḍiɗċ1=2µ#Ṫ

Try it online!

Uses 1-indexing (e.g. 1 -> 5, 2 -> 9, etc.), allowed by default on sequence challenges.

How it works

ŻḤḤ‘ḍḍiɗċ1=2µ#Ṫ - Main link. Takes n via STDIN
            µ#  - Execute the following on integers k = 1, 2, 3, ... until n integers return True:
Ż               -   Yield [0, 1, 2, ..., k]
 Ḥ              -   Unhalve; [0, 2, 4, ..., 2k]
  Ḥ             -   Unhalve; [0, 4, 8, ..., 4k]
   ‘            -   Increment; [1, 5, 9, ..., 4k+1]
       ɗ        -   Group the previous three links into a dyad.
                    Use l = [1, 5, ..., 4k+1] on the left and k on the right:
    ḍ           -     Each element in l is divisible by k?
      i         -     Index of k in l or 0?
     ḍ          -     1 or 0 is divisible by the index?
                    This yields a list with 2 1s for Hilbert primes
        ċ1      -   Count 1s
          =2    -   Equals 2?
              Ṫ - Take the last one i.e. the nth Hilbert prime

Husk, 8 bytes

!ü¦¡+4 5

Try it online! Uses 1-based indexing.

Explanation

!ü¦¡+4 5   Implicit input n.
   ¡       Iterate
    +4     addition of 4
       5   starting from 5: [5,9,13,17,..]
 ü         Uniquify by
  ¦        divisibility.
!          Get nth element.

Most of the work is done by ü¦. The function ü, when given a binary function f and a list, greedily constructs a subsequence where f x y is falsy for every (not necessarily adjacent) pair of elements x, y. In this case it picks a number if it's not divisible by any earlier pick. Then we just index into the resulting infinite list.

CJam, 36 33 32 23 bytes

5ri{_L+:L;{4+_Lf%0&}g}*

Try it online

The latest version is actually much more @MartinBüttner's than mine. The key idea in his suggested solution is to use two nested loops to find the n-th value that meets the condition. I thought I was being clever by using only a single loop in my original solution, but it turns out that the added logic cost more than I saved by not using a second loop.

Explanation

5       Push first Hilbert prime.
ri      Get input n and convert to integer.
{       Loop n times.
  _       Push a copy of current Hilbert prime.
  L       Push list of Hilbert primes found so far (L defaults to empty list).
  +       Prepend current Hilbert prime to list.
  :L      Store new list of Hilbert primes in variable L.
  ;       Pop list off stack.
  {       Start while loop for finding next Hilbert prime.
    4+      Add 4 to get next Hilbert number.
    _       Copy candidate Hilbert number.
    L       Push list of Hilbert primes found so far.
    f%      Element wise modulo of Hilbert number with smaller Hilbert primes.
    0&      Check for 0 in list of modulo values.
  }g      End while loop.
}*      End loop n times.

Julia, 73 bytes

n->(a=[x=5];while length(a)<n;x+=4;all(k->mod(x,k)>0,a)&&push!(a,x)end;x)

Thanks Alex A. for saving 11 bytes! This uses the same algorithm as the Matlab and Ruby answers. Since Julia arrays are one-indexed, this starts with f(1) == 5.

My first attempt, using the Lazy package, is 106 bytes. If you plan to run this in the REPL, make sure to add semicolons to the ends of the lines to suppress the infinite output. And call Pkg.Add("Lazy") if you don't already have it installed.

using Lazy
r=range
h=r(1,Inf,4)
p=@>>r() filter(n->n!=1&&all(map(x->mod(h[n],h[x])<1,2:n-1)))
f=n->h[p[n]]

Matlab, 74 83 bytes

function t=H(n)
x=5;t=x;while nnz(x)<n
t=t+4;x=[x t(1:+all(mod(t,x)))];end

Thanks to Tom Carpenter for removing 9 bytes!

Example use:

>> H(20)
ans =
   101

JavaScript (ES6), 73 bytes

n=>{for(i=0,t=2;i<=n;)i+=!/^(.(....)+)\1+$/.test(Array(t+=4));return t-1}

Just check Hilbert numbers one by one until we reach the nth Hilbert prime. Divisibility by Hilbert number is handled by regex.

Ruby, 60 bytes

h=->i{n=[];x=5;n.any?{|r|x%r<1}?x+=4: n<<x until e=n[i-1];e}

Only checks Hilbert prime factors.

Pyth, 21 bytes

Lh*4bye.fqZf!%yZyT1hQ

Try it online: Demonstration or Test Suite

Explanation:

Lh*4bye.fqZf!%yZyT1Q    implicit: Q = input number
L                       define a function y(b), which returns
 h*4b                      4*b + 1
                        this converts a index to its Hilbert number
       .f          hQ   find the first (Q+1) numbers Z >= 1, which satisfy:
           f      1        find the first number T >= 1, which satisfies:
            !%yZyT            y(Z) mod y(T) == 0
         qZ                test if the result is equal to Z 

                        this gives a list of indices of the first Q Hilbert Primes
      e                 take the last index
     y                  apply y and print

Mathematica, 65 bytes

Select[4Range[4^9]+1,Divisors[#][[2;;-2]]~Mod~4~FreeQ~1&][[#+1]]&

Generates the entire list and selects the element from it.

Minkolang 0.14, 46 37 32 bytes

I didn't realize that the gosub was totally unnecessary... >_>

n$z(xxi4*5+d(4-$d%)1=,z+$ziz-)N.

Try it here and check all test cases here.

Explanation

n$z                                 Take number from input and store it in the register
   (                                Open while loop
    xx                              Dump the stack
      i4*5+                         Loop counter times 4 plus 5 (Hilbert number)
           d                        Duplicate
            (                       Open while loop
             4-                     Subtract 4
               $d                   Duplicate stack
                 %                  Modulo
                  )                 Exit while loop when top of stack is 0
                   1=,              0 if 1, 1 otherwise
                      z             Push register value
                       +            Add
                        $z          Pop and store in register
                          iz-       Subtract z from loop counter
                             )      Exit while loop when top of stack is 0
                              N.    Output as number and stop.

The register is used to store the target index. The outer while loop calculates each Hilbert number and does some bookkeeping. The inner while loop checks each Hilbert number for primality. If a Hilbert number is not a Hilbert prime, then the target is incremented so that the outer while loop has to repeat (at least) one more time, effectively skipping Hilbert composites.

Haskell, 46 bytes

(foldr(\a b->a:[x|x<-b,mod x a>0])[][5,9..]!!)

An anonymous function.

The core is foldr(\a b->a:[x|x<-b,mod x a>0])[][5,9..], which iterates through the arithmetic progression 5,9,13,..., removing multiples of each one from the list to its right. This produces the infinite list of Hilbert primes. Then, !! takes the nth element.

I has tried making (\a b->a:[x|x<-b,mod x a>0]) pointfree but didn't find a shorter way.