Here's a way to find the number. It looks like it might not be the best (Sloane's A055315 has other formulas), but it works.
Let $n$ be the number of vertices and $e=n-1$ be the number of edges.
Claim: There is exactly one vertex of degree $3$, all others have degree $2$ or $1$. (Assume otherwise, and we violate the Handshaking Lemma, which stipulates the sum of degrees is $2e=2(n-1)$.)
So, the trees look like this:
If we delete the degree-$3$ vertex, we obtain three components. These three components only have vertices of degree $1$ and $2$, and so must be paths. Suppose the paths have $a$, $b$ and $c$ vertices. Then $a+b+c=n-1$ and $a,b,c$ are all positive integers. Conversely, given three paths of lengths $a \geq 1$, $b \geq 1$ and $c \geq 1$ and $a+b+c=n-1$, we can join their endpoints to a new vertex to create an $n$-vertex tree with exactly three $3$ leaves.
Let $a(n)$ be the number of partitions of $n-1$ into $3$ non-zero parts; this is equivalent to Sloane's A001399 which includes the formula $$a(n)=1+\left\lfloor\frac{(n-1)^2+6(n-1)}{12}\right\rfloor$$ among others.
Given an unlabeled tree, we can label the vertices in $n!$ (not necessarily distinct) ways. In this way, we generate $a(n)\ n!$ labeled trees.
In some cases (not many), there are symmetries which need to be accounted for. Because of symmetry, in this case
we count every labeling twice. This happens whenever $\{a,b,c\}$ have two parts of the same size.
And in this case
we count every labeling six times. This happens when $\{a,b,c\}$ has all three parts of the same size.
In these cases, we'll need to add a correction.
In the first case: we subtract $$\sum_{a \geq 1:n-1-a \text{ is even and } \geq 2} \tfrac{1}{2}n!=\frac{\lfloor (n-2)/2 \rfloor\ n!}{2}$$ since we have double-counted the $\frac{1}{2} n!$ graphs with $b=c=\tfrac{n-1-a}{2}$.
In the second case: if $n-1 \equiv 0 \pmod 3$, we first "uncorrect" the first correction, so we add back in $\tfrac{1}{2}n!$ and then subtract $\frac{5}{6}n!$ since we originally counted six times the $\frac{1}{6} n!$ graphs with $a=b=c$.
This gives $$\begin{cases} a(n)\ n!-\tfrac{\lfloor (n-2)/2 \rfloor\ n!}{2} & \text{if } n-1 \not\equiv 0 \pmod 3 \\ a(n)\ n!-\tfrac{\lfloor (n-2)/2 \rfloor\ n!}{2}+\tfrac{1}{2}n!-\tfrac{1}{6} n! & \text{if } n-1 \equiv 0 \pmod 3. \end{cases}$$ in agreement with Sloane's A055315.
Any graph on $n$ vertices with $k$ automorphisms – ways its vertices can be mapped onto themselves to preserve the edges – has $\frac{n!}k$ labellings. Applying this to the trees on seven vertices (the image below taken from Peter Steinbach's Field Guide to Simple Graphs, available on the OEIS under the links at A000055):
we see that the trees from left to right have the following numbers of automorphisms:
$$2,2,1,6,8,2,6,4,12,24,720$$
Dividing $7!=5040$ by each number gives the number of labellings of each of these trees:
$$2520,2520,5040,840,630,2520,840,1260,420,210,7$$
As expected, they add up to $7^5=16807$.
Finding the order of the automorphism group of a tree
As an example, take the second tree from the left. There is only one order-3 vertex, so it must stay fixed in any automorphism. Three paths radiate from this vertex – one of length 4 that must also stay fixed, and two of length 1 that can be swapped. Multiplying the numbers together gives two automorphisms.
For the rightmost star, while the central hub must stay fixed, the six spokes can be permuted in any order, yielding $6!=720$ automorphisms.
Best Answer
One approach is the following:
The result is going to give a uniformly random tree with the constraint, because in each $n$-vertex tree, there are $\binom{n+k}{n}$ ways to undo step 4, and $(n-1)n(n+1)(\dots)(n+k-2)$ ways to undo step 3 (for each vertex of degree 2, in order from smallest to largest, pick an edge to insert it on). So each $n$-vertex tree can be obtained in the same number of ways.
(I guess, equivalently, we could choose a random Prüfer code for an $(n+k)$-vertex tree, and delete the entries that only show up once, changing the labels of the rest in the same way as before. We'd have to be careful to distinguish labels that appeared $0$ times from labels that appeared $1$ time that we deleted, when we relabel.)
The tricky part is choosing $k$. Each label occurs exactly once in the Prüfer code with probability approximately $\frac 1e$, so the expected number of vertices of degree $2$ is approximately $\frac{n+k}{e}$. This makes it reasonable to choose $k = \lfloor \frac n{e-1}\rfloor$, for example.
The probability of getting exactly $k$ degree-$2$ vertices can be approximated by $\Pr[\text{Binomial}( n+k, \frac1e) = k]$, which is $O(\frac{1}{\sqrt n})$. This is not quite true, because the degrees are not independent, but still suggests that only $O(\sqrt n)$ trials are needed.
Here's a simple implementation using IGraph/M in Mathematica:
A faster, but somewhat less simple version: