Arturo, 34 bytes

f:$=>[x:&-1(x>0)?->+x*f x-1f x->1]

Try it!

The recurrence relation given in the question.

HOPS, 18 bytes

exp(x+x^2/2).*{n!}

Attempt This Online!

Japt, 16 12 11 bytes

1-indexed

È+T°*ZÔÅÎ}g

Try it

APL (Dyalog Unicode), 37 bytes

{⍵<2:1⋄+/(2⍵-1)×∇¨⍵-⍳2}

Try it online!

-1 byte from Adám.

Stax, 7 bytes

öƒ◙─v◙6

Run and debug it

K (ngn/k), 27 bytes

{*|x{x,+/(2#|x)*1,#2_x}/!2}

Try it online!

• x{...}/!2 setup a do-reduction with x (the input) as the number of iterations to run, seeded with 0 1 (the first two terms of the sequence)
• 1,#2_x build list containing how to multiply the previous values (i.e. T(n-1) by 1, and T(n-2) by n-1). Note that this infers the value of n from the length of the sequence so far.
• (2#|x)* get the last two values of the reversed sequence (i.e. (x[n-1];x[n-2]))
• x,+/ take their sum, and append to the sequence
• *| return the last value of the generated sequence

Husk, 7 bytes

#S=ηÖPḣ

Try it online! or Verify first 10 values

Same idea as the Jelly answer, and equally as slow.

0-indexed.

Explanation

#S=ηÖPḣ
      ḣ range 1..n
     P  permutations
#       number of arrays which satisfy:
 S      self
  =     equals
   η    indices of
    Ö   ordered elements?

k4, 53 bytes

{+/{*/[1+!x]*%*/[y#2]**/*/'1+(!x-2*y;!y)}[x]'!1+_x%2}

   {                                    }[x]'!1+_x%2 /pass inner lambda 2 args - x and each (') enumerate (!) 1+floor float-div 2
                             (!x-2*y;!y)             /enumerate
                           1+                        /add 1 to both list items
                        */'                          /multiply over each - this gives both denominator factorials
                      */                             /multiply them
              */[y#2]                                /make list of 2s of y-length and multiply over - gives y^2
                     *                               /multiply left and right args
             %                                       /reciprocal
    */[1+!x]                                         /x factorial
 +/                                                  /sum

examples:

  {+/{*/[1+!x]*%*/*/[y#2],*/'1+(!x-2*y;!y)}[x]'!1+_x%2}'[0 3 10 14]
1 4 9496 2390480f

Python 3, ^{81 \$\cdots\$ 57} 56 bytes

def f(n):
 a=b=i=1
 while i<n:a,b=a+i*b,a;i+=1
 return a

Try it online!

Uses 0 based indexing and implements the recursive formula.

Red, 53 bytes

f: func[n][either 1 > n: n - 1[1][n *(f n - 1)+ f n]]

Try it online!

Mathematica 46 42 bytes

f=2^(#/2)HypergeometricU[-#/2,.5,-.5]I^-#&

Try it online!

Thanks to mabel for helping me shave 4 bytes by using what I think is an anonymous function.

As seen in https://www.wolframalpha.com/input/?i=OEIS+A000085 and http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheSecondKind.html

Elm, 45 bytes

f n = if n<2 then 1 else (n-1)*f(n-2)+f(n-1)

Try it online!

Implementation of the recursive function.

JavaScript (ES6), 25 bytes

0-indexed. Uses the recurrence relation.

f=n=>n<2||f(--n)+n*f(n-1)

Try it online!

TI-BASIC, 28 22 bytes

∑((Ans nCr (2K))(2K)!/(2^KK!),K,0,int(.5Ans

Input is n in Ans.
Output is the nth telephone number.

I was going to use the first formula mentioned in the challenge, but testing it resulted in ERROR: OVERFLOWs even on smaller numbers because TI-BASIC handles n!! as two successive factorials instead of a double factorial.
At least it has a builtin for summation notation!

Update:

I found that \$(2k-1)!!=\frac{(2k)!}{2^kk!}\$, so i replaced the function i was using to \${{n}\choose{2k}}\frac{(2k)!}{2^kk!}\$, which is represented as (Ans nCr (2K))(2K)!/(2^KK!) in TI-BASIC.

Explanation:

∑((Ans nCr (2K))(2K)!/(2^KK!)               ;sum the equation mentioned above
                             ,K             ;using K as the loop variable
                               ,0           ;starting at 0
                                 ,int(.5Ans ;and ending at the floor of half of the input
                                            ;leave the result in Ans
                                            ;implicit print of Ans

Examples:

10:prgmCDGF26
            9496
5:prgmCDGF26
              26

Note: TI-BASIC is a tokenized language. Character count does not equal byte count.

GolfScript, 16 bytes

This is inspired by the Fibonacci GolfScript program.

~,1.@{)*1$+\}/\;

Try it online!

Explanation

~,               # Generate exclusive range from input to 0
  1.             # Make 2 copies of 1
    @            # Make the each target on the top
     {      }/   # Foreach over the range
      )          # Increment the current item
       *         # Multiply top by the current item
        1$+      # And add by the second-to-top
           \     # Swap items
              \; # Swap & discard the bottom item

Oasis, 7 6 5 bytes

àn*+1

Try it online!

# to compute the nth telephone number f(n):
à      # push the telephone numbers f(n-1) and f(n-2)
 n     # push n
  *    # multiply: n * f(n-2)
   +   # add: n * f(n-2) + f(n-1)

1      # base case: f(0) = 1

PHP, 48 bytes

function f($n){return$n<2?1:f(--$n)+$n*f($n-1);}

Try it online!

Implementation of the recursive function.

Original version: 54 bytes (port of @Noodle9 answer):

for($i=$j=1;++$n<$argn;$j=$k){$k=$i;$i+=$n*$j;}echo$i;

Try it online!

Japt, 28 26 bytes

_pZÌ+(ZÊÉ *ZgJÑ}g´U1õ ï)gJ

Try it

pZÌ+(ZÊÉ *ZgJÑ   make sequence with next element from previous sequence

_    ...  }g´U    repeat input -1 times 
  1õ ï)          starting with [[1,1]]
       gJ        implicit output last element of resulting sequence

C (gcc), 40 bytes

long f(long n){n=n<2?1:f(--n)+n*f(n-1);}

Try it online!

Merely copied from @Arnauld answer

Retina 0.8.2, 80 bytes

.+
$*_;;_;
{`^_(_*;)(_*;_+)
$1_$2,$2
+`(,_*)_(;_+;(_*))|,
$3$1$2
^;_*;(_*).*
$.1

Try it online! Explanation:

.+
$*_;;_;

Initialise the work area with n, i=0, and T=[1] (note that the list indices are reversed, so that the first element of T is T[i] and the last is T[0]) converted to unary.

{`

Loop n times (actually until the buffer stops changing, but see below).

^_(_*;)(_*;_+)
$1_$2,$2

Decrement n, increment i, and make extra copies of i and T[i].

+`

Repeat i times.

(,_*)_(;_+;(_*))|,
$3$1$2

Decrement the copy of i and add T[i-1] to the copy of T[i], stopping when the copy reaches zero by deleting the , entirely, causing the repeat to terminate. The result is therefore T[i+1]=T[i]+iT[i-1].

^;_*;(_*).*
$.1

When n is zero, replace the work area with T[i] converted to decimal. This causes the loop to terminate as the stages can no longer match.

J, 25 bytes

($:+]*$:@<:)@<:`1:@.(<&2)

Try it online!

$: is J's recursion operator, so this is a fairly straightforward impl of the recursive def.

I tried a couple other approaches (including using ^: power of operator) but wasn't able to get shorter than this.

Am curious if anyone can improve it...

Keg, 21 bytes

:1≤[_1|;:@Tƒ^:;@Tƒ*+]

Try it online!

Simply implements the formula in a recursive manner. The footer is mainly for testing purposes

Haskell, 29 bytes

f n|n<2=1|m<-n-1=f m+m*f(n-2)

Try it online!

Implements the recursive formula.

Python, 33 bytes

f=lambda n:n<2or~-n*f(n-2)+f(n-1)

Try it online!

Uses the recursive formula.

Jelly, 6 bytes

Œ!ỤƑ€S

A monadic Link accepting a non-negative integer which yields a positive integer.

Try it online! Or see the first 10 (n=10 is too slow)

How?

Builds all permutations of [1...n] (or if n=0 just [[]]) and counts the number which are involutions.

Œ!ỤƑ€S - Link: integer, n     e.g.  0        3
Œ!     - all permutations          [[]]     [[1,2,3],[1,3,2],[2,1,3],[2,3,1],[3,1,2],[3,2,1]]
    €  - for each:
   Ƒ   -   is invariant under:
  Ụ    -     grade                ( []       [1,2,3] [1,3,2] [2,1,3] [3,1,2] [2,3,1] [3,2,1] )
       -   }                        [1]     [1,      1,      1,      0,      0,      1]
     S - sum                        1       4

Charcoal, 20 bytes

⊞υ¹ＦＮ⊞υ⁺§υι×ι§υ⊖ιＩ⊟υ

Try it online! Link is to verbose version of code. Uses the recurrence relation. Explanation:

⊞υ¹

T(0)=1

ＦＮ

Loop n times.

⊞υ⁺§υι×ι§υ⊖ι

T(i+1)=T(i)+iT(i-1)

Ｉ⊟υ

Output T(n).

05AB1E, 12 8 bytes

-4 bytes by porting Stephen's cQuent answer

1λèsN<*+

Try it online!

Perl 6, 26 bytes

{{(1,1,++$×*+*...*)[$_]}}

Try it online!

cQuents, 10 bytes

=1:Z+Y($-2

1-indexed.

Try it online!

Explanation

=1           first term in sequence is 1
  :          given n, output nth term in sequence (1-indexed)
             each term is
   Z+                     (n-1) term +
     Y                                 (n-2) term *    
      ($-2                                          (index - 2)