How to generate heteroskedastic data for linear regression analysis given Y I have at m different points on a surface representing an organ n  measures of a organ property for n subjects (such as wall thickness). These values have been stored in a matrix Y with m columns and n rows. The measures at different points are highly correlated - the correlation coefficient between two columns of Y always ranges between 0.6 and 0.9. For each column $\tilde{Y}$ of Y I have computed a linear regression model of the form
$ \tilde{Y} = \beta X $
where X is a vector which contains n values of a variable such as age or height - one for each subject. X is always the same for all the m regression models. 
By doing this, I am trying to test at all the points of the surface where there is a significant association between the values at that point and the variable in X - age for example. By using this approach (mass univariate analysis) I could discover the regional effects of that clinical variable on the organ. However, as the number of points under study is greater than 100k, I have to apply to the p-values associated to each correlation coefficient a multiple testing correction that I have specifically created for this problem.
Unfortunately, the data are heteroscedastic and therefore one of the linear regression assumptions is violated and I would like to test how the failing of this assumptions affects the results that I have obtained. In particular, in order to test that I would like to generate a vector X with no relationship ($\beta=0$) with $ \tilde{Y}$ and that would make the variance of Y unequal along the range of X (heteroskedasticy). But in doing so, I would like to maitain untouched the values of Y as I believe that the correlation between its columns plays an important role.
Do you have any idea/suggestion on how I can generate such data, please?
I already succeed to generate heteroskedastic data by generating X using a normal distribution and by adding to each column of Y an additional term generated with a normal distribution with mean 0 an variance equal to X, but in this way I am losing the correlation between the columns of Y.
 A: There's a big gain to be gotten by reorganizing your data. Right now, you're treating the wall thickness of each point in the heart as a completely separate dependent variable (DV), when clearly there are meaningful relationships among all these things. Instead, move the position information to $X$, so now you have only one DV. (How to encode the position in $X$ is potentially a deep topic, for which see a textbook on spatial data analysis, but whatever coordinate information you use to identify the positions in the first place is a decent place to start.) Note that you now need to add a subject identifier to $X$ since you now have $m$ rows for each subject instead of just one. Now, you can build one big regression model instead of $m$ small ones. You can look at, e.g., the overall effect of age on wall thickness by using a main effect of age, or the effect of age on wall thickness at a particular point on the heart by using an interaction term. You will probably want to use a mixed model where each person gets a random intercept, and perhaps also each position gets a random intercept.
Heteroscedasticity may or may not show up again, which there is no doubt many ways to respond to, but you can avoid it optimistically biasing your conclusions by predictively validating your model instead of relying on significance testing.
A: +1 to @Kodiologist, that's clearly the answer as best we can tell from what has been written.  Moreover, this is by definition an XY problem, since how to generate heteroscedastic data is what you want to do to "test how the failing of this assumptions affects the results that I have obtained".  
That said, I can somewhat address the topic you have specified.  Namely, you cannot generate $X$ data to create a null relationship in which the existing $Y$ data will be heteroscedastic.  You will not be able to do this because heteroscedasticity is in the $Y$ data, and not in the $X$ data.  That's why you have been able to generate new $Y$ data with heteroscedasticity, but have been having trouble replicating the heteroscedasticity by generating pseudo-random $X$ data.  To understand this more fully, it may help you to read my answer here: What does having “constant variance” in a linear regression model mean?
A: Based on this post: sk learn quantile regression
In python you could do this as follows:
import numpy as np
import matplotlib.pyplot as plt

rng = np.random.RandomState(42)
x = np.linspace(start=0, stop=10, num=100)
X = x[:, np.newaxis]
y_true_mean = 10 + 0.5 * x
y_heteroskedastic = y_true_mean + rng.normal(loc=0, scale=0.5 + 0.5 * x, size=x.shape[0])

plt.scatter(X,y_heteroskedastic, label='heteroskedastic data')
plt.plot(X,y_true_mean, '-g', label='true mean')
plt.legend()
plt.show()


