The basic tool for this is lm; note that aov is a wrapper for lm.

In particular, if you have some grouping variable (factor), $g$, and a continuous covariate $x$, the model y ~ x + g would fit a main effects ANCOVA model, while y ~ x * g would fit a model which includes interaction with the covariate. aov will take the same formulas.

Pay particular attention to the Note in the help on aov.

As for + vs *, russellpierce pretty much covers it, but I'd recommend you look at ?lm and ?formula and most especially section 11.1 of the manual An Introduction to R that comes with R (or you can find it online if you haven't figured out how to find it on your computer; most easily, this involves finding the "Help" pull down menu in either R or RStudio).

I recommend getting and reading Discovering Statistics using R by Field. He has a nice section on ANCOVA.

To run ANCOVA in R load the following packages:

car
compute.es
effects
ggplot2
multcomp
pastecs
WRS

If you are using lm or aov (I use aov) make sure that you set the contrasts using the "contrasts" function before doing either aov or lm. R uses non-orthogonal contrasts by default which can mess everything up in an ANCOVA. If you want to set orthogonal contrasts use:

contrasts(dataname$factorvariable)=contr.poly(# of levels, i.e. 3) 

then run your model as

model.1=aov(dv~covariate+factorvariable, data=dataname)

To view the model use:

Anova(model.1, type="III") 

Make sure you use capital "A" Anova here and not anova. This will give results using type III SS.

summary.lm(model.1) will give another summary and includes the R-sq. output.

posth=glht(model.1, linfct=mcp(factorvariable="Tukey"))  ##gives the post-hoc Tukey analysis
summary(posth) ##shows the output in a nice format.

If you want to test for homogeneity of regression slopes you can also include an interaction term for the IV and covariate. That would be:

model=aov(dv~covariate+IV+covariate:IV, data=dataname)

If the interaction term is significant then you do not have homogeneity.

Here is a complementary documentation http://goo.gl/yxUZ1R of the procedure suggested by @Butorovich. In addition, my observation is that when the covariate is binary, using summary(lm.object) would give same IV estimate as generate by Anova(lm.object, type="III").

We use Regression analysis to create models which describe the effect of variation in predictor variables on the response variable. Sometimes if we have a categorical variable with values like Yes/No or Male/Female etc. the simple regression analysis give multiple results for each value of the categorical variable. In such scenario, we can study the effect of the categorical variable by using it along with the predictor variable and comparing the regression lines for each level of the categorical variable. Such an analysis is termed as Analysis of Covariance also called as ANCOVA.

Example
Consider the R built-in data set mtcars. In it we observer that the field am represents the type of transmission (auto or manual). It is a categorical variable with values 0 and 1. The miles per gallon value (mpg) of a car can also depend on it besides the value of horse power (hp). We study the effect of the value of am on the regression between mpg and hp. It is done by using the aov() function followed by the anova() function to compare the multiple regressions.

Input Data
Create a data frame containing the fields mpg, hp and am from the data set mtcars. Here we take mpg as the response variable, hp as the predictor variable and am as the categorical variable.

input <- mtcars[,c("am","mpg","hp")]
head(input)

When we execute above code, it produces following result:

                  am  mpg  hp
Mazda RX4          1 21.0 110
Mazda RX4 Wag      1 21.0 110
Datsun 710         1 22.8  93
Hornet 4 Drive     0 21.4 110
Hornet Sportabout  0 18.7 175
Valiant            0 18.1 105

ANCOVA Analysis
We create a regression model taking hp as the predictor variable and mpg as the response variable taking into account the interaction between am and hp.

Model with interaction between categorical variable and predictor variable

Create regression model1

result1 <- aov(mpg~hp*am,data=mtcars)
summary(result1)

When we execute above code, it produces following result:

            Df Sum Sq Mean Sq F value   Pr(>F)    
hp           1  678.4   678.4  77.391 1.50e-09 ***
am           1  202.2   202.2  23.072 4.75e-05 ***
hp:am        1    0.0     0.0   0.001    0.981    
Residuals   28  245.4     8.8                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

This result shows that both horse power and transmission type has significant effect on miles per gallon as the p-value in both cases is less than 0.05. But the interaction between these two variables is not significant as the p-value is more than 0.05.

Model without interaction between categorical variable and predictor variable

Create the regression model2

result2 <- aov(mpg~hp+am,data=mtcars)
summary(result2)

When we execute above code, it produces following result:

            Df Sum Sq Mean Sq F value   Pr(>F)    
hp           1  678.4   678.4   80.15 7.63e-10 ***
am           1  202.2   202.2   23.89 3.46e-05 ***
Residuals   29  245.4     8.5                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

This result shows that both horse power and transmission type has significant effect on miles per gallon as the p-value in both cases is less than 0.05.

Comparing two models
Now we can compare the two models to conclude if the interaction of the variables is truly insignificant. For this we use the anova() function.

 anova(result1,result2)

 Model 1: mpg ~ hp * am
 Model 2: mpg ~ hp + am
   Res.Df    RSS Df  Sum of Sq     F Pr(>F)
 1     28 245.43                           
 2     29 245.44 -1 -0.0052515 6e-04 0.9806

As the p-value is greater than 0.05 we conclude that the interaction between horse power and transmission type is not significant. So the mileage per gallon will depend in a similar manner on the horse power of the car in both auto and manual transmission mode.