How can I get slope and standard error at several levels of a continuous by continuous interaction in R? I'm comparing the slopes of several different response variables (DVs; representing different populations) to a set of predictors (IVs).  For some DVs a 2-way interaction (continuous by continuous) is supported.  To facilitate my comparison of IV coefficients I'd like to plot the slope estimates and 95% CI on a single graph (separate graph for each IV), and for the DV's with an interaction I'd like to plot the slope at ~3 values of the continuous moderator variable (e.g., "DV 1" in figure below).

I'm sure there is a variety of ways to get these values, but I'm hoping someone can point me to a simple bit of code or a package that can help automate this process for me.  I should also note my models are from lme4.
The 'effects' package handily calculates predicted values at user-specified levels of the moderator variable, but doesn't provide slopes or SE to my knowledge (although I could figure these out from the predicted values, I'm hoping for a more stream-lined method).  
Here is some toy data, although it doesn't produce an interaction like I show in the figure;
set.seed(50)
x1 <- rnorm(100,2,10)
x2 <- rnorm(100,2,10)
y1 <- x1+x2+x1*x2+rnorm(100,0,100)

model1<-lm(y1 ~ x1*x2)

And here is the predicted values plotted from 'effects', but I want the slopes and SE of these lines... 
library(effects)
model1.eff<-effect("x1*x2",model1,xlevels=3)
plot(model1.eff,multiline=T,ci.style="bands")
as.data.frame(model1.eff)

 A: While @wools answer appears more than adequate, here is another alternative that allows the calculation of marginal effect of x1 given x2, from a single model output without centering the x variables;
According to http://statistics.ats.ucla.edu/stat/r/faq/concon.htm ;
where the model is
y ~ β0 + β1x1 + β2x2+ β3x1∗x2
then slope for x1 at a given value of x2 is β1 + β3 * x2
So I can choose a few values of x2 as;
at.x2<-c(-6, 1, 6)

slopes <- coef(model1)["x1"] + coef(model1)["x1:x2"] * at.x2

According to How to calculate the standard error of the marginal effects in interactions (robust regression)?
Standard error for slopes = sqrt(var(b1) + var(b3) x2^2 + 2 x2 * cov(b1,b3) ) 
estvar<-vcov(model1); model1.vcov<-as.data.frame(as.matrix(estvar))
var.b1<-model1.vcov["x1","x1"]
var.b3<-model1.vcov["x1:x2","x1:x2"]
cov.b1.b3<-model1.vcov["x1","x1:x2"]

SEs <- rep(NA, length(at.x2))
for (i in 1:length(at.x2)){
  j <- at.x2[i]  
  SEs[i] <- sqrt(var.b1 + var.b3 * j^2 + 2*j* cov.b1.b3)
}

cbind(SEs, slopes)

