Why is this sequence of recurrence relevant?
Calculus of differences

We make Kaya's nice observations slightly more formal. The sequence $T^k$ is the $k$th difference sequence of the original sequence $T$. You can think of it as a discrete $k$th derivative. For example, if $T$ is a $k$th degree polynomial then $T^k$ is constant. We will see below what happens when $T$ is exponential.

For a sequence $T$ (all sequences range over the positive integers), let $\Delta(T)$ be defined by $$ \Delta(T)(n) = T(n+1)-T(n). $$
So $T^k = \Delta^{(k)}(T)$, that is $T^k$ is obtained by applying $k$ times the operator $\Delta$. The operator $\Delta$ has an inverse, which corresponds to indefinite integration:
$$ T(n) = T(1) + \sum_{k=1}^{n-1} \Delta(T)(k). $$
If we are only given the sequence $\Delta(T)$ then we don't know $T(1)$, and we can think of it as a constant of integration. For concreteness, define another operator $\Sigma$ given by
$$ \Sigma(T)(n) = \sum_{k=1}^{n-1} T(k). $$
The "fundamental theorem of calculus" then states that
$$ \Delta(\Sigma(T)) = T $$
and
$$ \Sigma(\Delta(T)) = T - T(1). $$
Another useful property of the operators $\Delta$ and $\Sigma$ is their linearity.

Furthermore, we have the nice formula $$\Delta^{(k)}(n^k) = k!,$$ which can be proved by induction on $k$. This formula also implies that $$\Delta^{(k+1)}(n^k) = 0.$$ In fact, this is true for any polynomial of degree at most $k$. Conversely, $\Sigma(n^k)$ is some polynomial of degree $k$.

We can do the same for exponentials:
$$ \Delta(c^n) = (c-1)c^n, \quad \Sigma(c^n) = \frac{c^n-1}{c-1}.$$

Finally, it will be useful to use the shift operator $S$ given by $$S(T)(n) = T(n+1).$$
This operator commutes with $\Delta$: $S(\Delta(T)) = \Delta(S(T))$.

Solving the recurrence

Now to the problem at hand. The sequence $T$ satisfies $$S(T) = 2T + (n+1)^2. $$
Taking the third difference, we get $$S(\Delta^{(3)}(T)) = 2\Delta^{(3)}(T).$$
Here we used the linearity of $\Delta$, the commutativity of $\Delta,S$, and the fact that the third derivative of $(n+1)^2$ vanishes. We can conclude that
$$ \Delta^{(3)}(T) = C 2^n, \quad C = \Delta^{(3)}(T)(1)/2. $$
We now want to repeatedly apply the operator $\Sigma$. Each time we apply $\Sigma$, we get some constant of integration. After the first time, we get
$$ \Delta^{(2)}(T) = C 2^n + C_1 $$
for some $C_1$, using $\Sigma(2^n) = 2^n-1$. While we can express $C_1$ in terms of $T$, it is better not to. The next time we apply $\Sigma$, $C_1$ turns into a linear terms, and so
$$ \Delta^{(1)}(T) = C 2^n + C_1 n + C_2 $$
for some $C_2$. The next time we apply $\Sigma$, $C_1 n$ turns into a quadratic term, and so
$$ T = C 2^n + P_2 $$
for some quadratic polynomial $P_2$.

It remains to find $C$ and $P_2$. Recall that $C = \Delta^{(3)}(T)(1)/2$. Kaya calculated $\Delta^{(3)}(T)(1) = 12$ and so $C = 6$. Therefore
$$ T = 6 \cdot 2^n + P_2. $$
The first few values in $T$ are $1,6,21$, and so $P_2(1) = -11$, $P_2(2) = -18$, $P_2(3) = -27$. Write $P_2 = B_2 n^2 + B_1 n + B_0$, and so $\Delta^{(2)}(P_2) = 2B_2$. On the other hand, $\Delta(-11,-18,-27,\ldots) = (-7,-9,\ldots)$ and so $\Delta^{(2)}(-11,-18,-27,\ldots) = (-2,\ldots)$, and we conclude that $B_2 = -1$. We can now calculate $B_1 n + B_0 = (-10,-14,\ldots)$, and so $B_1 = -4$ and $B_0 = -6$. We can conclude that $P_2 = -n^2 - 4n - 6$, and we get Kaya's formula.

Generalities

More generally, if $T$ is a sequence satisfying $T(n+1) = c T(n) + P_k(n)$ for some polynomial $P_k$ of degree $k$, then the very same method shows that
$$ T(n) = C c^n + P_{k+1}(n) $$
for some constant $C$ and polynomial $P_{k+1}$ of degree $k+1$. In order to find $C$ and $P_{k+1}$, we can either follow the route used above - express $C$ in terms of $c$, $k$ and $\Delta^{(k+1)}(T)(1)$ - or calculate $T(1),\ldots,T(k+2)$ and solve linear equations in the unknowns $C$ and the coefficients of $P_{k+1}$.

Even more generally, we can handle recurrences of the form $T(n) = c_1 T(n-1) + \cdots + c_l T(n-l) + P_k(n)$. As before, after applying $\Delta^{(k+1)}$, the inhomogeneous terms $P_k$ disappears, and we can solve the resulting recurrence relation. Applying $\Sigma^{(k+1)}$, we will get an "error term" of the form $P_{k+1}$ as before.

General recurrence relations have solutions which are sums of expressions of the form $n^r c^n$ for integer $r \geq 0$, so we also need to calculate the (iterated) discrete integrals of these functions; the result of applying $\Sigma$ once is $P_{r+1}(n) c^n$ for some polynomial $P_{r+1}$ of degree $r+1$, and so applying it $k+1$ time we obtain $P_{r+k+1}(n) c^n$ for some polynomial $P_{r+k+1}$ of degree $r+k+1$. As before, we can find the coefficients of all polynomials by solving linear equations, given enough values of the sequence $T$.