I think what you want is just:
fit <- lm(y ~ x+ f1*f2, data=dfrm)
This will give a different prediction for each level of the interaction of f1 with f2. It's just one model but could be "queried" for predictions with the predict
function using any desired combo of f1 and f2. You should look at ?formula and spend some time understanding how linear models get interpreted.