SAKO, 87 81 bytes

G=1
*1)DRUKUJ(9,0):G
G=(48×I*2+64×I+20)×G/(3×I*2+11×I+10)
POWTORZ:I=0(1)-1
KONIEC

SAKO isn't really a competitive language here, still was worth a try.
Outputs the infinite sequence with no input, by printing.
E notation will be used for numbers with length above 9 digits.

It's possible to lower the byte count to 74 with the code below, but floating-point inaccuracies show up as fast as G(8).

G=1
*1)DRUKUJ(9,0):G
G=(140/(3×I+5)-84/(I+2)+16)×G
POWTORZ:I=0(1)-1
KONIEC

Japt -N, 32 bytes

¶Tªß´U *[20I48]xV=È*UpYÃ/[AB3]xV

Try it

¶Tªß´U *[20I48]xV=È*UpYÃ/[AB3]xV     :Implicit input of integer U
¶                                    :Is equal to
 T                                   :0
  ª                                  :Logical OR with
   ß                                 :  Recursive call with argument
    ´U                               :    Prefix decrement U
       *                             :  Multiply by
        [20I48]                      :    Array of integers, where I=64
               x                     :    Reduce by addition after
                  È                  :    Mapping each X at 0-based index Y through the following function
                   *                 :      Multiply X by
                    UpY              :        U raised to the power of Y
                V=                   :    Assign that function to variable V
                       Ã             :    End reduce
                        /            :  Divide by
                         [AB3]       :    Array of integers, where A=10 & B=11
                              xV     :    Reduce by addition after mapping through V

The following computes Gessel walks of length 2n

import numpy as np

def walk(n) :

    if n == 1:

        return(['E','W','I','O'])

    else:

        A = []

        for w in walk(n-1):

            A = A  + [w+'E',w+'W',w+'O',w+'I']

        return(A)

def weight(s):

    if s == 'E':

        return([1,0])

    if s == 'W':

         return([-1,0])

    if s == 'O':    

         return([1,1])

    if s == 'I' :   

        return([-1,-1])

    else:

        d = [0,0]

        for w in s :

            d = np.add(d,weight(w))

    return(d)
    
def pwalk(w):

    a = [0,0]

    A = []

    for s in w:

        a = np.add(a,weight(s))

        if (a < [0,0]).any():

            return(False)

        else:

            continue

    return(True)       

 def Gessel(n):

    A = []

    for w in walk(2*n):

        if pwalk(w) and (weight(w)==[0,0]).all():

            A = A + [w]

    return(A)

The function Gessel(n) returns all the Gessel walks of length 2n.

Ruby, 46 bytes

f=->n{n<1?1:f[n-1]*(4-112*n+16*n*=3*n+5)/n+=2}

Try it online!

Charcoal, 30 bytes

→ＦＮＪ÷×ⅈ↨⊗⟦²⁴¦³²χ⟧ι↨⟦³¦¹¹χ⟧ι⁰Ｉⅈ

Try it online! Link is to verbose version of code. Outputs the \$n\$th value. Explanation: Uses a formula from the OEIS page.

→

Start with \$g(0)=1\$.

ＦＮ

Loop \$n\$ times.

Ｊ÷×ⅈ↨⊗⟦²⁴¦³²χ⟧ι↨⟦³¦¹¹χ⟧ι⁰

$$g(i+1)=\frac{48i^2+64i+20}{3i^2+11i+10}g(i)$$

Ｉⅈ

Output \$g(n)\$.

Actually counting the paths takes 70 bytes:

≔Ｅ²Ｅ²∧ιλυＦＮ«≔Ｅ⁺⟦Ｅυ⁰⟧υ⁺⟦⁰⟧κυＵＭυＥκ∧λ∧νΣ⟦§κ⊖ν§κ⊕ν§§υ⊖λ⊖ν§§υ⊕λ⊕ν⟧»ＩΣＥυΣＸι²

Try it online! Link is to verbose version of code. Explanation:

≔Ｅ²Ｅ²∧ιλυ

Start with 1 path from the origin to itself and 0 paths to 0, 1, 1, 0 and 1, 1, but with inverted indexing so that the origin is actually at the far corner of the lattice.

ＦＮ«

Loop n times.

≔Ｅ⁺⟦Ｅυ⁰⟧υ⁺⟦⁰⟧κυ

Add a protected boundary to the lattice so far. (This effectively allows Charcoal to index out of the lattice.) Note that this becomes part of the lattice for the next loop when the lattice grows a new boundary.

ＵＭυＥκ∧λ∧νΣ⟦§κ⊖ν§κ⊕ν§§υ⊖λ⊖ν§§υ⊕λ⊕ν⟧

For each square not on the boundary, sum the four neighbours with the allowed steps.

»ＩΣＥυΣＸι²

So far we've only calculated the number of routes to each square, but we also need to get back to the origin. Fortunately that's just the square of this, summed over all of the squares.

JavaScript (ES6), 48 bytes

Direct implementation of the recurrence formula given on OEIS.

f=n=>n?f(--n)*(20+64*n+48*n*n)/(10+11*n+3*n*n):1

Try it online!

JavaScript (ES6), 44 bytes

Derived from dharr's recursion formula, which is shorter.

f=n=>n?f(n-1)*(112*n-16*(n*=3*n+5)-4)/~-~n:1

Try it online!

JavaScript (ES6), 43 bytes

Also derived from dharr's formula. Doing the division first introduces small rounding errors.

f=n=>n?f(n-1)/(7*n+(n*=3*n-2)+2)*(16*n+4):1

Try it online!

CASIO BASIC (CASIO fx-9750GIII), 44 bytes

1
For 2→N To 84
2Ans(48N²-32N+4)÷(3N²+5N+2◢
Ans÷2
Next

prints until error (aka until the 84th term)

N>10 may be inaccurate because my calculator uses scientific notation for numbers above 9e9.

uses the recursive formula:

$$g(n)=\frac{g(n-1)(48n^2-32n+4)}{3n^2+5n+2}$$

Uiua 0.15.0-dev.2, ³⁰ 28 bytes ^SBCS

⧻⊚≡(/↧=0⊂⊃/↧⊣\+)⧅⋅⋅1:≡\+A₂×2

Try on Uiua Pad!

Takes \$n\$ and returns \$g(n)\$.

I thought there should be at least one solution that directly compute the paths.

Explanation

⧅⋅⋅1:≡\+A₂×2
          ×2 # 2n
     ≡\+A₂   # array of the four step vectors (0,1), (1,1), etc.
⧅⋅⋅1:        # choose 2n of them with replacement

⧻⊚≡(/↧=0⊃⊂/↧⊣\+)
   ≡(            ) # for each...
               \+  # running sum, coordinates after each step
          ⊂/↧⊣    # min of each column, and last row
         ⊃        # concatanate
     /↧=0          # is it all 0s?
⧻⊚               # how many TRUEs?

21 bytes

⁅/×/÷⊂16⊞+⊂⟜+₂10_3×6⇡

Try on Uiua Pad!

Direct computation using the formula in the question. ⁅ just rounds the final answer due to floating-point error.

APL(NARS), 38 chars

{0>r←⍵-1:1⋄(∇r)×16+(140÷5+3×r)-84÷2+r}

if one make the division on polynoms is possible reduce some chars because

÷/r⊥¨(48 64 20)(3 11 10) should be 16+(140÷5+3×r)-84÷2+r Test:

  {0>r←⍵-1:1⋄(∇r)×16+(140÷5+3×r)-84÷2+r}¨0,⍳10
1 2 11 85 782 8004 88044 1020162 12294260 152787976 1946310467 

this instead would be 32 chars lenght but fail the input 0

×/{16+(140÷5+3×⍵)-84÷2+⍵}¨0..⎕-1

Test:

×/{16+(140÷5+3×⍵)-84÷2+⍵}¨0..⎕-1
⎕:
  100x
24668399816955858855186657869247243539195399447597371872915574560359473983085405040263115861911922296227610827555120
  ×/{16+(140÷5+3×⍵)-84÷2+⍵}¨0..⎕-1
⎕:
  0
4
  ×/{16+(140÷5+3×⍵)-84÷2+⍵}¨0..⎕-1
⎕:
  2
11

R, 52 bytes

f=\(n)`if`(n,f(n-1)*(48*n^2-32*n+4)/(3*n^2+5*n+2),1)

Attempt This Online!

Using the recursion (see @dharr's answer).

R, 61 59 bytes

\(n,k=1:n-1)`if`(n,prod(16^n,5/6+k,.5+k)/prod(2+k,5/3+k),1)

Attempt This Online!

Gessel's formula.

Jelly, 19 bytes

6“¦€¤½‘÷+ⱮḶF×÷ƭƒ⁴*$

A full program that accepts a non-negative integer, \$n\$, and prints a positive float.

Inaccuracies creep in due to floating point arithmetic, and relatively quickly due to the order of operations (e.g. at \$n=15\$ the result is \$0.8\$ away). Still, it is theoretically accurate given enough precision.

Try it online! (The footer calls the program twice for each \$n\$ due to the Link itself not resetting the state of ƭ, which is why it's a full program.)

How?

6“¦€¤½‘÷+ⱮḶF×÷ƭƒ⁴*$ - Main Link: non-negative integer, N
 “¦€¤½‘             - [5, 12, 3, 10]
6      ÷            - divide by 6 -> [5/6, 2, 1/2, 5/3]
        +ⱮḶ         - add [0..N-1] -> [[5/6, 2, 1/2, 5/3], [11/6, 3, 3/2, 8/3], ...]
           F        - flatten -> P = [5/6, 2, 1/2, 5/3, 11/6, 3, 3/2, 8/3, ...]
                  $ - last two links as a monad - f(N):
                ⁴   -   16
                 *  -   exponentiate {N} -> 16^N
               ƒ    - starting with {16^N} reduce {P} by:
              ƭ     -   alternate between:
            ×       -     multiply; and
             ÷      -     divide

Integer calculation version, 21 bytes

Ḷ×6+€“¦€“¤½‘×€4,1PP:/

A monadic Link that accepts a non-negative integer and yields a positive integer.

Try it online!

How?

Ensure all multiplications are done first and a single integer division is performed at the end.

Ḷ×6+€“¦€“¤½‘×€4,1PP:/ - Link: non-negative integer, N
Ḷ                     - lowered range {N} -> [0..N-1]
 ×6                   - times six -> [0, 6, 12, ..., 6N-6]
   +€“¦€“¤½‘          - add [[5, 12], [3, 10]] to each
            ×€4,1     - multiply each by [4, 1]
                 P    - product
                  P   - product
                   :/ - reduce {this pair of integers} by integer division

Brute force version, 20 bytes

1pØ.;N$ṗḤ+\€»ƑƇ0ṪÐṂL

A monadic Link that accepts a non-negative integer and yields a positive integer.

How?

1pØ.;N$ṗḤ+\€»ƑƇ0ṪÐṂL - Link: non-negative integer, N
1pØ.                 - 1 Cartesian product [0,1] -> [[1, 0], [1, 1]]
    ;N$              - concatenate its negative
                         -> Deltas = [[1, 0], [1, 1], [-1, 0], [-1, -1]]
        Ḥ            - double {N} -> 2N
       ṗ             - {Deltas} Cartesian power {2N} -> AllWalkDeltas
         +\€         - cumulative reduce each by addition -> AllWalks
            »ƑƇ0     - keep those invariant under max with zero -> QuadrantRestictedWalks
                ṪÐṂL - count of those with minimal tails

MATL, 24 bytes

Thanks to Kevin Cruijssen for a correction in the second formula.

1i:"'!Q15'd@ZQ*[I5H]@ZQ/

Try at MATL online!

This inputs \$n\$ and outputs \$g(n)\$. The code is based on the recurrence relation given by OEIS,

$$g(n+1)=\frac{48n^2+64n+20}{3n^2+11n+10}g(n),$$

An equivalent expression, which is more convenient in this case due to MATL's \$1\$-based indexing, is

$$g(n)=\frac{48n^2-32n+4}{3n^2+5n+2}g(n-1).$$

How it works

1          % Push 1
i          % Input: n
:          % Range: gives [1 2 ... n]
"          % For each k in [1 2 ... n]
  '!Q15'   %   Push this string
  d        %   Consecutive differences (of code points). Gives [48 -32 4]
  @        %   Push k
  ZQ       %   Evaluate polynomial [48 -32 4] at k
  *        %   Multiply
  [I5H]    %   Push array [3 5 2]
  @        %   Push k
  ZQ       %   Evaluate polynomial [3 5 2] at k
  /        %   Divide
           % End (implicit)
           % Display stack (implicit)

05AB1E, 24 bytes

λ48Nn*32N*-4+*3Nn*5N*+Ì÷

Using the recursive formula:

$$g(n)=\frac{g(n-1)(48n^2-32n+4)}{3n^2+5n+2}$$ where \$g(0)=1\$.

Outputs the infinite sequence.^†

Try it online.

25 bytes

1ÝD!‚D(«I·ãεη€øO¤_s˜d«P}O

Using a hard-coded approach of actually calculating all Gessel walks.

Given \$n\$, outputs \$g(n)\$. Times out for \$n\geq5\$.

Try it online.

34 33 bytes

16Im5 6/Iݨ©+P2z®+PPI>!5 3/®+P*/ò

Using the given formula of the challenge description:

$$g(n)=\frac{16^n(5/6)_n(1/2)_n}{(2)_n(5/3)_n}$$ where \$(a)_n\$ is \$a(a+1)\cdots(a+n-1)\$.

Given \$n\$, outputs \$g(n)\$.^†

Try it online.

^{†: Results start to differ between the two for \$n\geq17\$, possibly due to floating point inaccuracies. 🤔}

Explanation:

λ               # Start a recursive environment,
                # to output the infinite sequence
                # Starting implicitly at g(0)=1
                # and where every following g(n) is calculated by:
 48Nn*32N*-4+   #  (48n²-32n+4)
             *  #  Multiplied by the implicit previous term g(n-1)
          ÷     #  Divided by:
 3Nn*5N*+Ì      #   (3n²+5n+2)

1Ý              # Push [0,1]
  D!            # Duplicate it, factorial on each value: [1,1]
    ‚           # Pair the two together: [[0,1],[1,1]]
     D(         # Duplicate it, negate each inner value: [[0,-1],[-1,-1]]
       «        # Merge the two together: [[0,1],[1,1],[0,-1],[-1,-1]]
        I·      # Push double the input
          ã     # Cartesian product to get all possible 2n-sized lists
ε               # Map over each:
 η              #  Pop and push its prefixes
  €ø            #  Map over each prefix: zip/transpose, swapping rows/columns
    O           #  Sum each inner-most list together
     ¤          #  Push the last pair (without popping the list)
      _         #  Check for both whether they're equal to 0
                #  (aka, they ended back at the origin)
     s          #  Swap so the list is at the top of the stack again
      ˜         #  Flatten it
       d        #  Check for each intermediate value that it's non-negative >=0
                #  (aka, we've stayed in the non-negative quadrant)
        «       #  Merge this list of checks to the earlier pair
         P      #  Check if all are truthy by taking its product
}O              # After the map: sum the amount of truthy result
                # (which is output implicitly as result)

16Im            # Push 16 to the power of the input
5 6/            # Push 5/6
    Iݨ         # Push a list in the range [0,input-1]
       ©        # Store this list in variable `®` (without popping)
        +       # Add 5/6 to each value in this list
         P      # Take the product
2z              # Push 1/2
  ®+P           # Calculate its rising factorial again using list `®`
P               # Take the product of all three values on the stack
I>!             # Push the factorial of input+1
                # (which equals the rising factorial of 2)
5 3/®+P*        # Multiply it to the rising factorial of 5/3
/               # Divide the two numbers from one another
 ò              # Round it to the nearest integer
                # (after which the result is output implicitly)

Nekomata + -n, 12 bytes

Äᵐ{įŋ∫}∫Ɔžᵚ≥

Attempt This Online!

Äᵐ{įŋ∫}∫Ɔžᵚ≥
Ä               Multiply the input by 2
 ᵐ{   }         Map over the list [0, 1, 2, ..., 2n-1]
   į                Push either [0,1] or [1,0]
    ŋ               Optionally negate the list
                        (Possible results: [0,1], [0,-1], [1,0], [-1,0])
     ∫              Cumulative sum
                        (Possible results: [0,1], [0,-1], [1,1], [-1,-1])
       ∫        Cumulative sum
        Ɔ       Split the list into the initial elements and the last element
         ž      Check if the last element is 0
          ᵚ≥    Check if the initial elements are all greater than or equal to 0

-n counts the number of solutions.

☾, 19 chars (38 bytes)

󰲡⅚!ˣ16ˣ½!ˣ÷2!ˣ⅗􍨅!ˣ▢

Note: ☾ uses a custom font, but you can run this code in the Web UI

Google Sheets, 88 bytes

=LET(a,LAMBDA(a,n,LET(x,n-1,IF(n=0,1,(20+64*x+48*x^2)*a(a,x)/(10+11*x+3*x^2)))),a(a,A1))

Maple, 51 bytes

f:=n->f(n-1)*(48*n^2-32*n+4)/(3*n^2+5*n+2);f(0):=1;

OEIS formula (10+11*n+3*n^2)*a(n+1) = (20+64*n+48*n^2)*a(n) modified to give a(n) from a(n-1).

APL+WIN, 73 bytes

Prompts for n

0⍕((2*((2÷3)+4×n))××/!((4 1 5÷3 2 6)+0,n,n)-1)÷(○,1)××/!((n←⎕)+(5÷3),2)-1

Try it online! Thanks to Dyalog Classic