Desmos, 57 bytes

f(1)=1
f(n)=max(1,∑_{d=2}^{n^½}0^{mod(n,d)}f(d)f(n/d))

Uses the given formula.

Try It On Desmos!

Try It On Desmos! - Prettified

It's really annoying that I have to define a base case f(1)=1 because otherwise Desmos would complain that I don't have a base case for my recursive function, even though f(1) is literally impossible to call recursively given how I defined f. In fact an explicit base case isn't even needed at all, since the base cases are already dealt with within the definition of f. I could literally put a completely absurd base case like f(-100000)=100000000000, for example, and all of the test cases would still work.

R, 70 bytes

f=\(x,y=2,z=0)`if`(x<y^2,max(z,1),f(x,y+1,z+`if`(x%%y,0,f(y)*f(x/y))))

Try it online!

Ugly, but well...

x=n is the real unknown of the function.

y is a potential divider of x. We increment it through recursive calls from 2 to ceiling(sqrt(x)).

z is an accumulator, counting the possible factorizations we have found with dividers smaller than y.

I wanted to golf it much more, with 2 variables instead of 3, but primes and squares of primes are annoying when you cheat too hard. Or maybe there's a smart trick I haven't found.

05AB1E, 29 20 17 15 bytes

<λèλNÑÂø2äн¦èPO

-9 bytes thanks to a tip of @emanresuA to use an alternative approach
-3 bytes thanks to @emanresuA directly

Try it online or verify the infinite list, starting at n=0.

Explanation:

 λ           # Start a recursive environment
< è          # To output the (implicit) (input-1)'th term
             # (which will be output implicitly afterwards)
             # (Staring implicitly with a(0)=1)
             # Where every following a(n) is calculated by:
   λ         #  Push the list of previous terms [a(0),a(1),...,a(n-1)]
    N        #  Push n
     Ñ       #  Pop and push a list of its divisors
      Â      #  Bifurcate it; short for Duplicate & Reverse copy
       ø     #  Create pairs with the values in the two lists
        2ä   #  Split it into two (almost†) equal-sized halves
          н  #  Pop and keep the first halve
           ¦ #  Remove the first [1,n] pair
    è        #  0-based index each inner value into the previous terms `λ`
     P       #  Take the product of each inner pair
             #  (if n is a prime number, the empty list will become 1)
      O      #  Them sum all of those together
             #  (for prime numbers, the 1 will remain 1)

† I say almost equal-sized halves, since square values will be odd-length lists (e.g. the pairs of 16 are [[1,16],[2,8],[4,4],[8,2],[16,1]], resulting in 'halves' [[[1,16],[2,8],[4,4]], [[8,2],[16,1]]]).

Google Sheets, 83 bytes

lambda(f,x,y,z,if(x<y*y,max(z,1),f(f,x,y+1,z+if(mod(x,y),,f(f,y,2,)*f(f,x/y,2,)))))

Anonymous lambda function that uses recursion. Based on Evargalo's R answer.

Call the function with the function itself as the first argument, \$n\$ as second argument, \$2\$ third, blank or zero fourth.

Ungolfed:

=let( 
  f, lambda(self, x, y, z, 
    if(x < y * y, 
      max(z, 1), 
      self(self, x, y + 1, z + if(mod(x, y), 
        0, 
        self(self, y, 2, 0) * self(self, x / y, 2, 0) 
      )) 
    ) 
  ), 
  f(f, A1, 2, 0) 
)

Ruby, 53 54 56 bytes

f=->n{[(2..n).sum{|x|x>n/x||n%x>0?0:f[x]*f[n/x]},1].max}

Try it online!

Previous attempt based on a false assumption, now fixed.

Wolfram Language (Mathematica), 77 71 bytes

Try it online!

a@n_?PrimeQ:=1
a@n_:=Sum[a[d]*a[n/d],{d,Select[Divisors@n,2<=#<=n/#&]}]

SageMath, 91 87 bytes

a=lambda n:int(is_prime(n))or sum(a(d)*a(n//d)for d in range(2,n) if n%d==0 and d*d<=n)

Python apart from the function 'is_prime'.

Shortened with Neil's advice.

Albert.Lang's hint removes 'is_prime':

Python, 74 bytes

a=lambda n:max(1,sum(a(d)*a(n//d)for d in range(2,n)if d*d<=n and n%d==0))

Haskell, 98 #96 92 bytes

s=z#1
z=0:z
(h:t)#i|n<-max 1h=n:zipWith(+)t(do p<-tail$take i s++z;(0<$[2..i])++[p*n])#(i+1)

Attempt This Online!

-2 remembering that thing I ALWAYS forget about

-4 thanks to xnor LMAO

Outputs infinitely. Goofy DP/sieve kind of thing, adding factorization counts up cumulatively.

s=z#1
z=0:z

z is an infinite list of zeroes (i.e. 0 prepended to itself). The sequence s is the function (#) applied to z and 1.

(h:t)#i|n<-max 1h

(#) destructures its first argument into its first element h and remainder t, and its second argument is called i. n is h, unless h is 0, in which case it's 1.

=n: [...] #(i+1)

(#) returns a list where the first value is n, followed by itself called on an incremented i and...

zipWith(+)t( [...] )

every element of t added to the corresponding element of...

do p<-tail$take i s++z;

the concatenation of, for every p in the first i elements of s except the first (infinitely padded with z, so zipWith doesn't truncate)...

(0<$[2..i])++[p*n]

i - 1 zeroes followed by the product of p and n.

Vyxal, 13 bytes

λKḢṪvxḂ*Ih∑1∴

Try it Online! Kinda ugly but it works.

λ             # Define a recursive function f taking n
    vx        # Do a recursive call on each of 
 K            # the divisors of n
  ḢṪ          # aside from 1 and n
      Ḃ*      # Multiply each f(k) by f(n/k)
        Ih    # Take the first half
          ∑   # and sum it
           1∴ # or return 1 if it's empty

Haskell, 49 bytes

f n=max 1$sum[f a*f b|a<-[2..n],b<-[2..a],a*b==n]

Try it online!

Jelly, 14 bytes

ḊŒċP⁼¥Ƈ⁸߀€ZPS

A monadic Link that accepts a positive integer and yields the count of factorisation trees.

Try it online!

How?

ḊŒċP⁼¥Ƈ⁸߀€ZPS - Link: positive integer, N
Ḋ              - dequeue -> [2..N]  (an empty list for N=1)
 Œċ            - unordered pairs with replacement -> [[2,2],[2,3],...,[2,N],[3,3],...,[N,N]]
                                                     (an empty list for N=1)
      Ƈ        - keep those Pairs for which:
     ¥ ⁸       -   last two links as a dyad - f(Pair, N):
   P           -     product
    ⁼          -     equals {N}?
        ߀€    - call this Link for each value of each remaining pair
           Z   - transpose  (an empty list for N=1)
            P  - product (vectorises)  (product of an empty list is 1)
             S - sum

JavaScript (Node.js), 57 54 50 bytes

f=(n,i=1)=>n%1?0:n>=++i*i&&f(n,i)+f(n/i)*f(i)||i<3

Try it online!

If there's no way to factor, then it's a prime and return 1

Charcoal, 35 bytes

Ｆ⊕Ｎ⊞υ∨ΣＥΦ…²ι¬∨›×κκι﹪ικקυκ§υ÷ικ¹Ｉ⊟υ

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

Ｆ⊕Ｎ

Loop n+1 times.

⊞υ∨ΣＥΦ…²ι¬∨›×κκι﹪ικקυκ§υ÷ικ¹

Calculate the next a(n) from its nontrivial minor factors, unless there are none, in which case assume a(n)=1. Note that Sum of an empty list is actually None in Charcoal, but it doesn't matter here because that's still falsy.

Ｉ⊟υ

Output the nth entry. (If you wanted all the entries from 0 to n then you could remove the ⊟.)

If you're willing to accept floating-point inaccuracy, then for 31 bytes:

Ｆ⊕Ｎ⊞υ∨ΣＥΦ…·²₂ι¬﹪ικקυκ§υ÷ικ¹Ｉ⊟υ

Try it online! Link is to verbose version of code. Explanation: Iterates up to the floating-point square root instead of filtering on the square being not greater than the number being factored.

In practice the floating-point library can take square roots accurately enough for values up to 2¹⁰⁶ which is far more than Charcoal will be able to iterate over before the heat death of the universe.