Reports

As @NelsonGon mentioned, vcov() works. Please see the below example using Swiss data.

data(swiss)

### multiple linear model, swiss data
lmod <- lm(Fertility ~ ., data = swiss)

vcov(lmod)

The covariance matrix as below:

                 (Intercept)   Agriculture   Examination    Education      Catholic Infant.Mortality
(Intercept)      114.6192408 -0.4849476484 -1.2025734658 -0.281265331 -0.0221836036    -3.2658448131
Agriculture       -0.4849476  0.0049426416  0.0043708713  0.004789532 -0.0005112844     0.0065656539
Examination       -1.2025735  0.0043708713  0.0644541409 -0.027310637  0.0051339487     0.0003482484
Education         -0.2812653  0.0047895318 -0.0273106371  0.033499469 -0.0030003666     0.0122667258
Catholic          -0.0221836 -0.0005112844  0.0051339487 -0.003000367  0.0012431162    -0.0027467320
Infant.Mortality  -3.2658448  0.0065656539  0.0003482484  0.012266726 -0.0027467320     0.1457098919

79497352