I would recommend using the new M12 function `AsymptoticSolve`

for this. Your equation:

```
eqn = FD1((d-2)/2, ηs) + FD1((d-2)/2, ηs - vd) == 2 FD1((d-2)/2, η0);
```

We need to find the zero order approximation of `ηs`

when `vd`

it is small:

```
Simplify(Solve(eqn /. vd -> 0), (η0 | vd) ∈ Reals)
```

{{ηs -> η0}}

Now use `AsymptoticSolve`

:

```
AsymptoticSolve(eqn, {ηs, η0}, {vd, 0, 5})
```

{{ηs ->

vd / 2 + ((-6 + d) (-2 + d) (-1 + d) vd ^ 4) / (

1536 η0 ^ 3) + ((2 – d) vd ^ 2) / (16 η0) + η0}}

**Appendix**

If you have an earlier version of Mathematica, so you do not have access to `AsymptoticSolve`

, you can try using the cloud instead. For example, define:

```
asymptoticSolve(args__) := CloudEvaluate(System`AsymptoticSolve(args))
```

Then use `asymptoticSolve`

instead of `AsymptoticSolve`

.