Denote $s_n=(q^n)_\infty=\prod_{k\geqslant n}(1-q^k)$. Note that $$\sum_{k\geqslant n}q^ks_{k+1}=\sum_{k\geqslant n}(s_{k+1}-s_k)=1-s_n.\quad\quad(\clubsuit)$$ The idea is to reduce your sum to such telescoping sums. For $k=0,1,\ldots,n$ denote $$ c_k:=(-1)^{k+1}{n\choose k}_qq^{k(k+1)/2}(q^{n+1})_{n-k}, $$ where ${n\choose k}_q=\frac{(q)_n}{(q)_k(q)_{n-k}}$ is a $q$-binomial coefficient. I claim that $$ (-1)^{n+1}q^{(n+1)(3n+2)/2}x(qx)_n\\=c_0q^{2n+1}x+c_1q^{2n}x(q^{2n+1}x)_1+c_2q^{2n-1}x(q^{2n}x)_2+\ldots+c_{n}q^{n+1}x(q^{n+2}x)_n. \quad \quad (\heartsuit)$$

First of all, I show how $(\heartsuit)$ yields your formula. Applying $(\heartsuit)$ for $x=q^j$ and using $(\clubsuit)$, your sum reads as $$ S=\sum_{i=0}^{n}c_i(1-s_{2n+1-i}). $$ I claim that $$-\sum_{i=0}^n c_is_{2n+1-i}=s_1=\sum_{k\in \mathbb{Z}} (-1)^kq^{k(3k-1)/2}\quad\quad(\spadesuit)$$ (the last equality is Euler's Pentagonal theorem). Thus $S-s_1=\sum_{i=0}^n c_i$, which is easy to see to be a polynomial in $q$ of degree $n(3n+1)/2$. On the other hand, $S$ is divisible by $q^{(n+1)(3n+2)/2}$ (that is seen from its very definition), thus this polynomial must be equal to $-s_1 \pmod {q^{(n+1)(3n+2)/2}}=-\sum_{|k|\leqslant n} (-1)^kq^{k(3k-1)/2}$, and we get your formula for $S$.

It remains to prove $(\heartsuit)$ and $(\spadesuit)$. I start with $(\spadesuit)$. Since $s_{2n+1-i}(q^{n+1})_{n-i}=s_{n+1}$, $(\spadesuit)$ reads as $$ s_{n+1}\sum_{i=0}^n (-1)^{i+1}{n\choose i}_q q^{i(i+1)/2}=-s_1, $$ or, equivalently, $$ \sum_{i=0}^n (-1)^{i}{n\choose i}_q q^{i(i+1)/2}=(q)_{n}, $$ which is a partial case of the $q$-binomial theorem $(x)_n=\sum_{i=0}^n (-1)^{i}{n\choose i}_q q^{i(i-1)/2}x^i$ for $x=q$.

Well, it remains to prove $(\heartsuit)$. For this we rewrite $(q^{2n+2-k}x)_k$ as $(-1)^kx^kq^{k(4n+3-k)/2}(q^{-2n-1}x^{-1})_k$ and get some variant of $q$-Vandermonde convolution identity.

Let ${n\choose k}_q=\frac{(q)_n}{(q)_k(q)_{n-k}}$ denote a $q$-binomial coefficient. We start with the following version of $q$-Vandermonde convolution identity: $$ (x-y)(x-qy)\ldots(x-q^{n-1}y)\\=\sum_{k=0}^n(-1)^{k}q^{k(k-1)/2}{n\choose k}_q(y)_{n-k}(xq^{1-k})_k\quad\quad(\diamondsuit) $$ A short proof of $(\diamondsuit)$ may be obtained by checking it when $x=q^i$, $y=q^{-j}$ for non-negative integers $i,j$ such that $i+j\leqslant n$. In these points the values of both sides of $(\diamondsuit)$ are equal to 0, unless $i+j=n$ and we take the summand with $k=i$ in the right hand side. At this point, the identity is staightforward. It remains to observe that this collection of points is enough for reconstruction of the polynomial of degree at most $n$.

We apply $(\diamondsuit)$ for $y=q^{n+1}$, $x=tq^{2n+1}$ and multiply both sides by $-q^{2n+1}$. This gives, denoting for $k=0,1,\ldots,n$ $$ c_k:=(-1)^{k+1}{n\choose k}_qq^{k(k+1)/2}(q^{n+1})_{n-k}, $$ the following identity: $$ (-1)^{n+1}q^{(n+1)(3n+2)/2}(qt)_n\\=c_0q^{2n+1}+c_1q^{2n}(q^{2n+1}t)_1+c_2q^{2n-1}(q^{2n}t)_2+\ldots+c_{n}q^{n+1}(q^{n+2}t)_n. \quad \quad (\heartsuit)$$

Denote now $s_n=(q^n)_\infty=\prod_{k\geqslant n}(1-q^k)$. Note that $$\sum_{k\geqslant n}q^ks_{k+1}=\sum_{k\geqslant n}(s_{k+1}-s_k)=1-s_n.\quad\quad(\clubsuit)$$ The idea is to reduce your sum to such telescoping sums.

Applying $(\heartsuit)$ for $t=q^j$ and multiplying by $q^j$ and using $(\clubsuit)$, your sum reads as $$ S=\sum_{i=0}^{n}c_i(1-s_{2n+1-i}). $$ I claim that $$-\sum_{i=0}^n c_is_{2n+1-i}=s_1=\sum_{k\in \mathbb{Z}} (-1)^kq^{k(3k-1)/2}\quad\quad(\spadesuit)$$ (the last equality is Euler's Pentagonal theorem). Thus $S-s_1=\sum_{i=0}^n c_i$, which is easy to see to be a polynomial in $q$ of degree $n(3n+1)/2$. On the other hand, $S$ is divisible by $q^{(n+1)(3n+2)/2}$ (that is seen from its very definition), thus this polynomial must be equal to $-s_1 \pmod {q^{(n+1)(3n+2)/2}}=-\sum_{|k|\leqslant n} (-1)^kq^{k(3k-1)/2}$, and we get your formula for $S$.

It remains to prove $(\spadesuit)$. Since $s_{2n+1-i}(q^{n+1})_{n-i}=s_{n+1}$, $(\spadesuit)$ reads as $$ s_{n+1}\sum_{i=0}^n (-1)^{i+1}{n\choose i}_q q^{i(i+1)/2}=-s_1, $$ or, equivalently, $$ \sum_{i=0}^n (-1)^{i}{n\choose i}_q q^{i(i+1)/2}=(q)_{n}, $$ which is a partial case of the $q$-binomial theorem $(x)_n=\sum_{i=0}^n (-1)^{i}{n\choose i}_q q^{i(i-1)/2}x^i$ for $x=q$.