BUG: quantile with list of quantiles fails on MultiIndex column and groupby
- [x] I have checked that this issue has not already been reported.

- [x] I have confirmed this bug exists on the latest version of pandas.

- [ ] (optional) I have confirmed this bug exists on the master branch of pandas.

---

**Note**: Please read [this guide](https://matthewrocklin.com/blog/work/2018/02/28/minimal-bug-reports) detailing how to provide the necessary information for us to reproduce your bug.

#### Code Sample, a copy-pastable example

```python
import pandas as pd
import numpy as np

steps = 6
simulations = 20
alpha = 0.05

idx_interim_date_forecast = pd.IndexSlice["2020-04-01":]

df_return_forecast = pd.concat(
    [
        pd.DataFrame(
            data=np.random.rand(steps, simulations),
            index=pd.date_range(
                start=idx_interim_date_forecast.start, periods=steps, freq="M"
            ),
            columns=pd.MultiIndex.from_product(
                [("sample",), range(0, simulations)], names=["scenario", "simulation"]
            ),
        ),
        pd.DataFrame(
            data=np.random.rand(steps, simulations),
            index=pd.date_range(
                start=idx_interim_date_forecast.start, periods=steps, freq="M"
            ),
            columns=pd.MultiIndex.from_product(
                [("trend",), range(0, simulations)], names=["scenario", "simulation"]
            ),
        ),
    ],
    axis=1,
    sort=True,
)

df_return_forecast.groupby(axis=1, level=0).quantile(q=alpha / 2)
df_return_forecast.groupby(axis=1, level=0).quantile(q=1 - alpha / 2)
df_return_forecast.groupby(axis=1, level=0).quantile(q=[alpha / 2, 1 - alpha / 2])

```
```python
Traceback (most recent call last):

  File "<ipython-input-1-6a99a78036cb>", line 37, in <module>
    df_return_forecast.groupby(axis=1, level=0).quantile(q=[alpha / 2, 1 - alpha / 2])

  File "C:\Users\Kurt\Anaconda3\lib\site-packages\pandas\core\groupby\groupby.py", line 1951, in quantile
    indices = np.arange(len(result)).reshape([len(q), self.ngroups]).T.flatten()

```

#### Problem description

I would expect the ```quantile``` call with an iterable list of quantiles to return at the specified locations. Individual calls at the required locations returns the correct quantiles. The issue seems to be a result of the multiindex column/groupby.

If instead I do the following I am able to get the desired result without error. 
```python
df_return_forecast.T.groupby(axis=0, level=0).quantile(q=[alpha / 2, 1 - alpha / 2]).T

``` 
So there is a workaround but the behaviour is not as expected and I am not sure of the performance issues it my induce for large dataframes.

#### Expected Output
As per result obtained by:
```python
df_return_forecast.T.groupby(axis=0, level=0).quantile(q=[alpha / 2, 1 - alpha / 2]).T

```
#### Output of ``pd.show_versions()``

<details>

INSTALLED VERSIONS
------------------
commit           : None
python           : 3.7.6.final.0
python-bits      : 64
OS               : Windows
OS-release       : 10
machine          : AMD64
processor        : Intel64 Family 6 Model 60 Stepping 3, GenuineIntel
byteorder        : little
LC_ALL           : None
LANG             : en
LOCALE           : None.None

pandas           : 1.0.3
numpy            : 1.18.1
pytz             : 2019.3
dateutil         : 2.8.1
pip              : 20.0.2
setuptools       : 46.1.3.post20200325
Cython           : 0.29.16
pytest           : 5.4.1
hypothesis       : 5.10.4
sphinx           : 3.0.2
blosc            : None
feather          : None
xlsxwriter       : 1.2.8
lxml.etree       : 4.5.0
html5lib         : 1.0.1
pymysql          : None
psycopg2         : None
jinja2           : 2.11.2
IPython          : 7.13.0
pandas_datareader: None
bs4              : 4.9.0
bottleneck       : 1.3.2
fastparquet      : None
gcsfs            : None
lxml.etree       : 4.5.0
matplotlib       : 3.2.1
numexpr          : 2.7.1
odfpy            : None
openpyxl         : 3.0.3
pandas_gbq       : None
pyarrow          : 0.16.0
pytables         : None
pytest           : 5.4.1
pyxlsb           : None
s3fs             : None
scipy            : 1.4.1
sqlalchemy       : 1.3.15
tables           : 3.6.1
tabulate         : None
xarray           : None
xlrd             : 1.2.0
xlwt             : 1.3.0
xlsxwriter       : 1.2.8
numba            : 0.48.0

</details>
