We propose simultaneous mean-variance regression for the linear estimation and approximation of conditional mean functions. In the presence of heteroskedasticity of unknown form, our method accounts for varying dispersion in the regression outcome across the support of conditioning variables by using weights that are jointly determined with the mean regression parameters. Simultaneity generates outcome predictions that are guaranteed to improve over ordinary least-squares prediction error, with corresponding parameter standard errors that are automatically valid. Under shape misspecification of the conditional mean and variance functions, we establish existence and uniqueness of the resulting approximations and characterize their formal interpretation and robustness properties. In particular, we show that the corresponding mean-variance regression location-scale model weakly dominates the ordinary least-squares location model under a Kullback-Leibler measure of divergence, with strict improvement in the presence of heteroskedasticity. The simultaneous mean-variance regression loss function is globally convex and the corresponding estimator is easy to implement. We establish its consistency and asymptotic normality under misspecification, provide robust inference methods, and present numerical simulations that show large improvements over ordinary and weighted least-squares in terms of estimation and inference in finite samples. We further illustrate our method with two empirical applications to the estimation of the relationship between economic prosperity in 1500 and today, and demand for gasoline in the United States.