Using @BrianMScott's recurrence, let's try to derive a formula for ourselves using linear algebra.
We have $\left[\begin{matrix}a_0 \\ a_1 \\ a_2\end{matrix}\right]=\left[\begin{matrix}1 \\ 1 \\ 2\end{matrix}\right]$. However, we want to make the beginning terms as simple as possible to make things easy for us. If you think about it, since each term is the sum of the last three, $a_{-1}$ must be $0$ in order for $a_2=a_1+a_0+a_{-1}$. Therefore, instead, we'll start with $\left[\begin{matrix}a_{-1} \\ a_0 \\ a_1\end{matrix}\right]=\left[\begin{matrix}0 \\ 1 \\ 1\end{matrix}\right]$
We want to figure out a matrix that gets us the next term by putting $a_1$ and $a_2$ at the top of the vector and then putting $a_3$, the new term, at the bottom. In order to put $a_1$ at the top, we make the first row $0 \ 1 \ 0$ so that the second term of the old vector will be the first term of the new vector. Similarly, to put $a_2$ in the middle, we make the second row $0 \ 0 \ 1$ so that the third term of the old vector will be the first term of the new vector. Finally, to make the new term the sum of the old terms, we make the last row $1 \ 1 \ 1$. Thus, we get the matrix:
$$\left[\begin{matrix}0 \ 0 \ 1 \\ 0 \ 1 \ 0 \\ 1 \ 1 \ 1\end{matrix}\right]$$
Each time we multiply the beginning vector by this matrix, we get a new term, so we have:
$$\left[\begin{matrix}0 \ 0 \ 1 \\ 0 \ 1 \ 0 \\ 1 \ 1 \ 1\end{matrix}\right]^n\left[\begin{matrix}0 \\ 1 \\ 1\end{matrix}\right]=\left[\begin{matrix}a_{n-1} \\ a_n \\ a_{n+1}\end{matrix}\right]$$
Now, we kind of have a formula for $a_n$ as the second term of the vector from this matrix-vector product, but we're taking a matrix to a power which is still really complicated. Therefore, in order to simplify this matrix-power, we need to diagonalize the matrix.
According to Wolfram Alpha, we have the characteristic polynomial $-\lambda^3+\lambda^2+\lambda+1$. Let's say the real root of this polynomial is $\alpha$. Thus, by synthetic division, we have:
$$-\lambda^3+\lambda^2+\lambda+1=(x-a)(-x^2+(1-a)x+(-a^2+a+1))$$
If we graph $-\lambda^3+\lambda^2+\lambda+1$ and then we graph the discriminant of this quadratic we have here, we see that the characteristic polynomial meets the x-axis when the discriminant is negative, meaning the other two roots of this polynomial are complex conjugates. This tells us that there are three distinct roots, one real and two complex, which we will call $\alpha$, $\beta$, and $\gamma$.
Since there are three distinct roots of the characteristic polynomial, we know the matrix is diagonizable. Now, we need to find the eigenvectors. Now, we know for any root $\lambda$, $-\lambda^3+\lambda^2+\lambda+1=0$, or that $1+\lambda+\lambda^2=\lambda^3$. Therefore, if we input the vector $\left[\begin{matrix}1 \\ \lambda \\ \lambda^2\end{matrix}\right]$, we get back $\left[\begin{matrix}\lambda \\ \lambda^2 \\ 1+\lambda+\lambda^2=\lambda^3\end{matrix}\right]$, making it an eigenvector of the eigenvalue $\lambda$. Thus, if we put all of our eigenvectors in a matrix, we get:
$$M=\left[\begin{matrix}1 \ 1 \ 1 \\ \alpha \ \beta \ \gamma \\ \alpha^2 \ \beta^2 \ \gamma^2 \end{matrix}\right]$$
Using Wolfram Alpha, we can find $M^{-1}$ in terms of $\alpha$, $\beta$, and $\gamma$:
$$M^{-1}=\left[\begin{matrix}\frac{\beta \gamma}{(\alpha-\beta) (\alpha-\gamma)} \ -\frac{\beta+\gamma}{(\alpha-\beta) (\alpha-\gamma)} \ \frac{1}{(\alpha-\beta) (\alpha-\gamma)} \\
\frac{\alpha \gamma}{(\beta-\alpha) (\beta-\gamma)} \ -\frac{\alpha+\gamma}{(\beta-\alpha) (\beta-\gamma)} \ \frac{1}{(\beta-\alpha) (\beta-\gamma)} \\
\frac{\alpha \beta}{(\gamma-\alpha) (\gamma-\beta)} \ -\frac{\alpha+\beta}{(\gamma-\beta) (\gamma-\beta)} \ \frac{1}{(\gamma-\alpha) (\gamma-\beta)}\end{matrix}\right]$$
With all of that information, we now have:
$$\left[\begin{matrix}0 \ 0 \ 1 \\ 0 \ 1 \ 0 \\ 1 \ 1 \ 1\end{matrix}\right]=M\left[\begin{matrix}\alpha \ 0 \ 0 \\ 0 \ \beta \ 0 \\ 0 \ 0 \ \gamma\end{matrix}\right]M^{-1}$$
and thus:
$$\left[\begin{matrix}a_{n-1} \\ a_n \\ a_{n+1}\end{matrix}\right]=\left[\begin{matrix}0 \ 0 \ 1 \\ 0 \ 1 \ 0 \\ 1 \ 1 \ 1\end{matrix}\right]^n\left[\begin{matrix}0 \\ 1 \\ 1\end{matrix}\right]=M\left[\begin{matrix}\alpha^n \ 0 \ 0 \\ 0 \ \beta^n \ 0 \\ 0 \ 0 \ \gamma^n\end{matrix}\right]M^{-1}\left[\begin{matrix}0 \\ 1 \\ 1\end{matrix}\right]$$
Compute the first matrix-vector product with $M^{-1}$:
$$\left[\begin{matrix}a_{n-1} \\ a_n \\ a_{n+1}\end{matrix}\right]=\left[\begin{matrix}0 \ 0 \ 1 \\ 0 \ 1 \ 0 \\ 1 \ 1 \ 1\end{matrix}\right]^n\left[\begin{matrix}0 \\ 1 \\ 1\end{matrix}\right]=M\left[\begin{matrix}\alpha^n \ 0 \ 0 \\ 0 \ \beta^n \ 0 \\ 0 \ 0 \ \gamma^n\end{matrix}\right]\left[\begin{matrix}\frac{1-\beta-\gamma}{(\alpha-\beta)(\alpha-\gamma)} \\ \frac{1-\alpha-\gamma}{(\beta-\alpha)(\beta-\gamma)} \\ \frac{1-\alpha-\beta}{(\gamma-\alpha)(\gamma-\beta)}\end{matrix}\right]$$
Now, if you use Wolfram Alpha to approximate the roots and then plug in $1-\alpha-\beta$ to a calculator, you'll find that it equals the other root. What this means is that $1-\alpha-\beta=\gamma$, $1-\alpha-\gamma=\beta$, and $1-\alpha-\beta=\gamma$. Thus, substitute:
$$\left[\begin{matrix}a_{n-1} \\ a_n \\ a_{n+1}\end{matrix}\right]=\left[\begin{matrix}0 \ 0 \ 1 \\ 0 \ 1 \ 0 \\ 1 \ 1 \ 1\end{matrix}\right]^n\left[\begin{matrix}0 \\ 1 \\ 1\end{matrix}\right]=M\left[\begin{matrix}\alpha^n \ 0 \ 0 \\ 0 \ \beta^n \ 0 \\ 0 \ 0 \ \gamma^n\end{matrix}\right]\left[\begin{matrix}\frac{\alpha}{(\alpha-\beta)(\alpha-\gamma)} \\ \frac{\beta}{(\beta-\alpha)(\beta-\gamma)} \\ \frac{\gamma}{(\gamma-\alpha)(\gamma-\beta)}\end{matrix}\right]$$
Do the next matrix-vector product with the diagonal matrix:
$$\left[\begin{matrix}a_{n-1} \\ a_n \\ a_{n+1}\end{matrix}\right]=\left[\begin{matrix}0 \ 0 \ 1 \\ 0 \ 1 \ 0 \\ 1 \ 1 \ 1\end{matrix}\right]^n\left[\begin{matrix}0 \\ 1 \\ 1\end{matrix}\right]=M\left[\begin{matrix}\frac{\alpha^{n+1}}{(\alpha-\beta)(\alpha-\gamma)} \\ \frac{\beta^{n+1}}{(\beta-\alpha)(\beta-\gamma)} \\ \frac{\gamma^{n+1}}{(\gamma-\alpha)(\gamma-\beta)}\end{matrix}\right]$$
Now, we want to find $a_n$, which is the second term, so instead of doing the whole matrix-vector product, we just need to multiple the second row of $M$, $\alpha \ \beta \ \gamma$ by this vector to get the formula:
$$\frac{\alpha^{n+2}}{(\alpha-\beta)(\alpha-\gamma)}+\frac{\beta^{n+2}}{(\beta-\alpha)(\beta-\gamma)}+\frac{\gamma^{n+2}}{(\gamma-\alpha)(\gamma-\beta)}$$
Note that this is very similar to one of the formulas found here, but because we said $a_0=a_1=1$ whereas they said $a_1=a_2=1$, we have the exponent $n+2$ where as the have the exponent $n+1$. Now, using a fairly accurate complex number arithmetic calculator (i.e. Wolfram Alpha) and the knowledge that the answer is an integer, we can calculate $a_n$ for any $n$ we want pretty quickly by plugging in the approximate values of $\alpha, \beta, \gamma$ and the value of $n$ into the formula. Hooray!
[1, 1, 2, 4, 7, 13, 24, 44, 81, 149, 274]
– Noble Mushtak May 08 '16 at 22:17