[R] heteroscedastic bivariate distribution with linear regression - prediction interval

Dear forum, which is the most suitable method to get the prediction interval of a bivariate normal distribution which is consistent with a linear model y = ax + b? I assume it is gls + predict. Am I correct? I'm rather new to R. Is there some reliable sample code for that problem? Thank you

[R] Single line command for variance of marginal distribution in 2 D linear regression?

Dear forum, we have a 2 dimensional normal distribution of (x, y) which is consistent with a linear regression modell of type y = ax + b. It has a heteroscedastic variance. Is there a straight and simple way with R to get the variance of the marginal distribution W(y|x=xi) where xi may be