Retina 0.8.2, 95 bytes

.+
$*a_b
{%`d
c$'¶$`cad
%`c
b$'¶$`cad
%`b
a$'¶$`ad$'¶$`da$'¶$`ada
}G`^(a)+_(?<-1>.)+$
\b(a+)_\1

Try it online! Link includes faster test cases. Explanation: Based on the Nitrodon's generator in his comment on the question.

.+
$*a_b

Convert the input to unary, insert a separator, and append the starting symbol (I'm using letters as some of the symbols have a special meaning.)

%`d
c$'¶$`cad

d can be substituted with either c or cad. (These should be (c) but I don't actually need the ()s.)

%`c
b$'¶$`cad

c can be substituted with either b or cad. (Again, omitting the ()s around the c.)

%`b
a$'¶$`ad$'¶$`da$'¶$`ada

b can be substituted with any of a, ad, da or ada. (The % limits the effect of $' and $` to a single line, so only the line with the match gets replicated.)

G`^(a)+_(?<-1>.)+$

Discard any substitutions which have become too long. This allows the loop to converge.

{`
}`

Repeat the above until no new substitutions can be made.

\b(a+)_\1

Count the number of exact matches.

HOPS, 23 bytes

A=1+A^2*x*(x+1);A+A^2*x

Attempt This Online!

Derived from the formula in Jonathan Allan's Jelly answer. Here A^2 is the generating function of S in Jonathan Allan's answer, which is OEIS A052705.

HOPS, 25 bytes

A=1+(1-A+(1+A+x*A)^2)*x/2

Attempt This Online!

Derived from the generating function formula on OEIS.

Wolfram Language (Mathematica), 68 64 bytes

SeriesCoefficient[(y=4t+4t^2)/(2-t-2t^2+(t+2)√(1-y)),{t,0,#}]&

Try it online! Uses the generating function given in OEIS, like dharr's answer, but SeriesCoefficient does the coefficient extraction for us. -4 bytes thanks to alephalpha who found an algebraically equivalent generating function that's shorter to use.

I realized that I can avoid the scientific-notation imperfection of the initial submission by replacing (...)^.5 with the equivalent √(...), since √ is a three-byte character for Mathematica.

Previous version (68 bytes):

SeriesCoefficient[(2-t-2t^2-(t+2)(1-4t-4t^2)^.5)/2/(t+1)^2,{t,0,#}]&

Nekomata + -n, 21 bytes

ŧĉᵐzç0ɔYṀƵ¿ᵗ{qCᵈƆM→>N

Attempt This Online!

A Rota-Baxter word with \$n\$ as can be represented as a list of \$n\$ non-negative integers, where each integer represents the level of nesting of the corresponding a. For example, a(a)a(a(a)a) can be represented as [0,1,0,1,2,1]. It can be easily verified that two different Rota-Baxter words have different representations.

Furthermore, the list representation should satisfy the following conditions:

1. Either the first or the last element is 0. So that the word either starts with a or ends with a.
2. There are no two consecutive elements that are equal. So that aa and )( does not occur.
3. For any contiguous subsequences of the list, say \$a_i, a_{i+1}, \ldots, a_{j-1}, a_j\$, we must have \$\max(a_i, a_j) + 1 \ge \min(a_{i+1}, \ldots, a_{j-1})\$. So that the word does not contain a balanced substring that is enclosed twice.

ŧĉᵐzç0ɔYṀƵ¿ᵗ{qCᵈƆM→>N
ŧ                       Find an n-tuple of 0, 1, 2, ..., n-1
 ĉᵐz                    such that no two consecutive elements are equal.
    ç0ɔ                 Prepend and append 0s.
       YṀƵ¿             Check that the list now contains consecutive equal elements,
                            which means that either the first or the last element
                            of the original list is 0.
           ᵗ{           Check that the following condition cannot be satisfied:
             q              Find a contiguous subsequence
              CᵈƆM          such that the max of the first and last elements
                  →         plus 1
                   >N       is less than all the other elements.

-n counts the number of solutions.

Maple,  72  66 bytes

n->(D@@n)(t->(2-t-2*t^2-(t+2)*sqrt(1-4*t-4*t^2))/2/(t+1)^2)(0)/n!;

The expression in t is the generating function given in OEIS. Differentiating n times, evaluating at t = 0 and dividing by n! gives the appropriate coefficient in the series.

Edit: shorter version above by using D for differentiation. Earlier:

n->eval(diff((2-t-2*t^2-(t+2)*sqrt(1-4*t-4*t^2))/2/(t+1)^2,t$n)/n!,t=0);

JavaScript (ES6), 83 bytes

-2 thanks to emanresu A

This is a port of Albert.Lang's Python answer, with even more recursion.

f=n=>n>1?2*h(n)+h(n-1):n
g=(n,k=n)=>--k>1?h(k)*g(n-k)+g(n,k):v=f(n)
h=n=>2*g(n-1)-v

Try it online!

R, 106 bytes

\(n,f=expression((2-F-2*F^2-(F+2)*(1-4*F-4*F^2)^.5)/2/(F+1)^2)){for(i in 1:n)f=D(f,"F")
eval(f)/prod(1:n)}

Attempt This Online!

Uses the generating function (2-t-2*t^2-(t+2)*sqrt(1-4*t-4*t^2))/(2*(t+1)^2) from the OEIS (slightly rearranged).
Hugely inspired by dharr's answer which led me to the R D=derivative function: upvote that.

R, 93 bytes

\(n)`if`(n<3,n,sum(2*?n,?n-1))
`?`=\(n,C=choose)sapply(1:(n/2),\(i,m=n-i)C(n,i-1)*C(2*m,n)/m)

Attempt This Online!

A port of the formula in Jonathan Allan's answer, which I haven't even tried to understand: upvote that too.

Charcoal, 50 bytes

⊞υ⁰≔⟦⁰⟧θＦＮ«≔∨¬ιζη≔↨υ⁰ζ≔Σ×θ⮌υε⊞θ⁺⁺εη⊗ζ⊞υ⁺↨θ⁰ε»Ｉ⁺η⊗ζ

Attempt This Online! Link is to verbose version of code. Explanation: Port of @Albert.Lang's ungolfed code, converted to dynamic programming for the functions e and j (b, o and n are only needed transiently.)

⊞υ⁰≔⟦⁰⟧θ

Start with e[0]=0 and j[0]=0.

ＦＮ«

Loop n times.

≔∨¬ιζη

b(m)=1 if m=1, otherwise it equals o(m-1), whose value would be still be in the variable z by this point on the second loop.

≔↨υ⁰ζ

o(m)=j[m-1]. (AtIndex(u, i) would also work here, but it wouldn't work below.)

≔Σ×θ⮌υε

n(m) is the dot product of e[:m] and j[:m][::-1].

⊞θ⁺⁺εη⊗ζ

e[m]=n(m)+b(m)+2*o(m).

⊞υ⁺↨θ⁰ε

j[m]=e[m]+n(m).

»Ｉ⁺η⊗ζ

Calculate the final number of Rota-Baxter words.

A port of a mixture of @xnor's port and @JonathanAllan's Jelly answer is only 42 bytes:

Ｉ∨ΣΣＥ⁻ＮＥ³¬ιＥ⊘ι÷Π⁻⊗⁻ι⊕λ…⁰ι×Π⊕⁺…⁰λ…⁰⁻ιλ⁻ι⊕λ¹

Attempt This Online! Link is to verbose version of code. Explanation: Charcoal has neither comb nor factorial built in, so I've broken the combinations down into products of ranges, except that the two n!s cancel, and instead of dividing two ranges from 1 to two values I can just take the product between those values.

Python 3.8, 98 bytes

Xnor's port of my Jelly answer.

lambda N:sum(comb(2*k,n)/k*comb(n,n+~k)for n in[N,N,N-1]for k in range(1,n))or 1
from math import*

An unnamed function that accepts a positive integer and returns the number of Rota-Baxter words with n as.

Try it online!

Note that for \$n>1\$ this returns a float, at \$n=28\$ this becomes inaccurate. An integer version that does not suffer from this is 100 bytes TIO:

lambda N:sum(comb(2*k,n)*comb(n,n+~k)//k for n in[N,N,N-1]for k in range(1,n))or 1
from math import*

Jelly,  28 27  19 bytes

_~⁹c×Ḥc¥:⁸ðþṫ-§ṚḄȯ1

A monadic Link that accepts a positive integer, \$n\$, and yields a positive integer, the number of Rota-Baxter words with \$n\$ as.

Try it online!

How?

The number of Rota-Baxter words with \$n\$ as with an a at only the left (or, equivalently, right) end is:

$$S(n) = \sum_{i=1}^{\lfloor \frac{n}{2} \rfloor}{\frac{\binom{n}{i-1}\binom{2(n-i)}{n}}{n-i}}$$

\$S(n-1)\$ also gives the number of Rota-Baxter words with \$n\$ as with an a at both ends, except when \$n=1\$ we have \$S(0)=0\$ rather than the required result of \$1\$.

The number of Rota-Baxter words with \$n>1\$ as is then \$2S(n)+S(n-1)\$.

_~⁹c×Ḥc¥:⁸ðþṫ-§ṚḄȯ1 - Link: integer, N
          ðþ        - {v in [1..N]} by {n in [1..N]} table of f(v,n):
_                   -   {v} subtract {n} -> v-n
 ~                  -   bitwise NOT  {that} -> n-v-1
  ⁹c                -   n choose {that} -> X = C(n, n-v-1)
       ¥            -   last two links as a dyad - f(v, n):
     Ḥ              -     double {v} -> 2v
      c             -     {that} choose {n} -> Y = C(2v, n)
    ×               -   {X} multiply {Y} -> C(n, n-v-1) × C(2v, n)
        :⁸          -   integer divide {that} by v
                                         -> C(n, n-v-1) × C(2v, n) ÷ v
            ṫ-      - tail from index -1 -> rows for n=N-1 and n=N
              §     - sums       -> [CountDoubleEnded, CountSingleEnded]
               Ṛ    - reverse    -> [CountSingleEnded, CountDoubleEnded]
                Ḅ   - from binary -> 2 × CountSingleEnded + CountDoubleEnded
                 ȯ1 - logical OR with 1 (replaces the 0 calculated for N=1 with 1)

Python, 116 bytes

f=lambda m:m>1and 2*j(m-1)+j(m-2)or m
n=lambda m:f(m)+sum(j(k)*n(m-1-k)for k in range(1,m-1))
j=lambda m:2*n(m)-f(m)

Attempt This Online!

How?

Here is an ungolfed version with an easier to understand recursion:

Python, 704 bytes

lambda m:2*o(m)+b(m)

# helper counts; all count balanced strings containing m a's subject to condition:
# e: enclosable: anything that may be inside a pair of parens within a valid RB
# b: bothsided:valid rb that starts and ends with an a (the same if m=1)
# n: nonesided: any enclosable string that has no a at either end; these are exactly the ones enclosed in a pair of non matching parens
# o: onesided: valid RB that has an a at exactly one end
# j: juxtaposable: anything balanced that may occur next to an a in a valid RB

e=lambda m:m>0 and n(m)+b(m)+2*o(m)
b=lambda m:m>0 and m<2 or j(m-2)
n=lambda m:m>2 and sum(e(k)*(j(m-1-k)) for k in range(m))
o=lambda m:m>1 and j(m-1)
j=lambda m:e(m)+n(m)

Attempt This Online!