I don't know how you performed your search, but there are a number of SAS procedures which employ generalized least squares estimation methods. Since you don't tell us anything about the problem other than that you want some way to deal with heteroskedasticity, it is a bit difficult to advise on what methods to employ. However, assuming that you are fitting a regression model in which the residual error is porportional to a power of the mean, then I would suggest looking at the MIXED procedure. Such estimation is discussed in the documentation of the REPEATED statement.
Actually, you might find it easier to use the NLMIXED procedure to estimate a model in which the residual variance is proportional to a power of the mean. NLMIXED code that fits a model in which the variance is proportional to the mean can be fit with code like the following:
proc nlmixed data=mydata;
mu = b0 + b1*x1 + b2*x2 + ...;
model y ~ normal(mu, mu*exp(2*log_sig));
estimate "Var(eps) for mu=1" exp(2*log_sig);
run;
There will be problems with the above model if your parameters produce a value of mu which is non-positive (zero or negative). A generalization of the above model which is consistent with the discussion of the POM model in the PROC MIXED documentation would be to fit the model:
proc nlmixed data=mydata;
mu = b0 + b1*x1 + b2*x2 + ...;
model y ~ normal(mu, (mu**theta)*exp(2*log_sig));
estimate "Var(eps) for mu=1" exp(2*log_sig);
run;
I made the assumption above that you were interested in a model in which the heteroskedasticity produced a residual variance which was proportional to a power of the mean. However, you might be interested in a model in which the residual variance differs according to the level of some categorical predictor variable. Such models can be handled quite easily with the MIXED procedure. Again, look at the documentation of the REPEATED statement and pay particular attention to the GROUP= option.
HTH,
Dale