I tried the following to check the identity for some specific $\lambda$ value, say
$$
\lambda = [4,3,2]\ .
$$
Then the L.H.S. in the expression above,
$$
s_{\lambda}(x_1 +y , x_2 +y ,\ldots , x_n +y)
$$
is computed by the following code:

```
R.<y> = PolynomialRing(QQ)
Sym = SymmetricFunctions(R)
Sym.inject_shorthands()
lam = Partition([3, 2, 1])
slam6 = s[lam].expand(6)
slam6.parent().inject_variables() # defining x0, x1, ... , x5
pol = slam6(*(xj + y for xj in (x0, x1, x2, x3, x4, x5)))
# pol # decomment to see it
```

Of course, we want to have it as a Schur function, taken from the above polynomial, yes, it is the following (manually split) expression:

```
sage: s.from_polynomial(pol)
896*y^6*s[] + 896*y^5*s[1] + 448*y^4*s[1, 1] + 112*y^3*s[1, 1, 1] \
+ 320*y^4*s[2] + 192*y^3*s[2, 1] + 48*y^2*s[2, 1, 1] + 32*y^2*s[2, 2] \
+ 8*y*s[2, 2, 1] + 40*y^3*s[3] + 24*y^2*s[3, 1] + 6*y*s[3, 1, 1] \
+ 4*y*s[3, 2] + s[3, 2, 1]
```

Let us check the coefficient in $\mu=[3,1]$ from the R.H.S. of the claimed identity.
We expect $24\;y^2\;s[3, 1]$ from the above expansion.

```
mu = Partition([3, 1])
skew = lam / mu
print(f'lam = {lam} , mu = {mu} , s(lam / mu) is {s(skew)}')
```

This gives

```
lam = [3, 2, 1] , mu = [3, 1] , s(lam / mu) is s[1, 1] + s[2]
```

We have to plug in now $y$ in the above Schur polynomial.
If my plug in uses the right conventions, then we obtain a different result.

```
sage: s(skew).expand(6)(y,y,y,y,y,y)
36*y^2
sage: s(skew).expand(1)(y)
y^2
```

Here are the Schur polynomials for the skew partitions $\lambda/\mu$ for further partitions $\mu$:

```
for k in [0..6]:
for mu in Partitions(k):
if lam.contains(mu):
print('lam = %s :: mu = %-9s :: s(lam / mu) = %s' % (lam, mu, s(lam/mu)))
```

This gives:

```
lam = [3, 2, 1] :: mu = [] :: s(lam / mu) = s[3, 2, 1]
lam = [3, 2, 1] :: mu = [1] :: s(lam / mu) = s[2, 2, 1] + s[3, 1, 1] + s[3, 2]
lam = [3, 2, 1] :: mu = [2] :: s(lam / mu) = s[2, 1, 1] + s[2, 2] + s[3, 1]
lam = [3, 2, 1] :: mu = [1, 1] :: s(lam / mu) = s[2, 1, 1] + s[2, 2] + s[3, 1]
lam = [3, 2, 1] :: mu = [3] :: s(lam / mu) = s[2, 1]
lam = [3, 2, 1] :: mu = [2, 1] :: s(lam / mu) = s[1, 1, 1] + 2*s[2, 1] + s[3]
lam = [3, 2, 1] :: mu = [1, 1, 1] :: s(lam / mu) = s[2, 1]
lam = [3, 2, 1] :: mu = [3, 1] :: s(lam / mu) = s[1, 1] + s[2]
lam = [3, 2, 1] :: mu = [2, 2] :: s(lam / mu) = s[1, 1] + s[2]
lam = [3, 2, 1] :: mu = [2, 1, 1] :: s(lam / mu) = s[1, 1] + s[2]
lam = [3, 2, 1] :: mu = [3, 2] :: s(lam / mu) = s[1]
lam = [3, 2, 1] :: mu = [3, 1, 1] :: s(lam / mu) = s[1]
lam = [3, 2, 1] :: mu = [2, 2, 1] :: s(lam / mu) = s[1]
lam = [3, 2, 1] :: mu = [3, 2, 1] :: s(lam / mu) = s[]
```