Showing that in a two variable model with an interaction, the lines for regression of $Y$ on $X_1$ at fixed values of $X_2$ cross at a point I am trying to show that in the model $E(Y) = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_1X_2$, the lines for regression of $Y$ on $X_1$ at various fixed values of $X_2$ all cross at a point. I've tried the following:
# Generate some data
x1 <- 11:30
x2 <- runif(20,5,95)

# True population coefficients
b0 <- 17
b1 <- 0.5
b2 <- 0.037
b3 <- -5.2
sigma <- 1.4

eps <- rnorm(x1,0,sigma)  # Error term
y <- b0 + b1*x1  + b2*x2  + b3*x1*x2 + eps  # Generate "true" y for our data

regression = lm(y ~ x1*x2)  # Regress y on x1, x2, and their interaction
coeffs = regression$coefficients  # Get their coefficients

temp = vector('list', 8)  # We want to plot 8 lines
for (X2 in 1:8) {  # Fixes values of X2
  points = coeffs[[1]] + coeffs[[2]]*1:8 + coeffs[[3]]*X2 + coeffs[[4]]*1:8*X2  # Y values for the lines we will plot
  temp[[X2]] = points
}

plot(1:8, temp[[1]], type='l')
lines(1:8, temp[[2]], type='l')

Clearly, though, these lines do not intersect, and so I am doing something wrong. Any help would be much appreciated!
 A: You have to assume $\beta_3\ne 0.$  When that's the case, the model can be rewritten as
$$\eqalign{y &= \beta_0 + \beta_1x_1 + \beta_2 x_2 + \beta_3 x_1x_2 \\
&= (\beta_0 - \beta_1\beta_2/\beta_3) + (\beta_1+\beta_3x_2)(x_1 + \beta_2/\beta_3) \\
&= \alpha + f(x_2)(x_1-\lambda)
}$$
where $\alpha=\beta_0 - \beta_1\beta_2/\beta_3,$ $\lambda = -\beta_2/\beta_3,$ and $f(x_2) = \beta_1 + \beta_3 x_2.$  These are all linear functions of $x_1$ and they all have the common value $\alpha$ for $x_1 = \lambda,$ showing they all intersect at the point $$(\lambda,\alpha) = (-\beta_2/\beta_3,\beta_0-\beta_1\beta_3/\beta_3).$$

When $\beta_3=0,$ the lines all intersect at the "point at infinity" given by the slope $\beta_1.$

Here is an illustration.  $\beta=(\beta_0, \ldots, \beta_3)$ is given in the title.  The lines correspond to various values of $x_2.$

If you would like to experiment, here's the R code to produce figures like it.
# set.seed(17)
beta <- rnorm(4)                        # Create random coefficients
lambda <- -beta[3]/beta[4]
alpha <- beta[1] + beta[2] * lambda
#
# Set up a plotting area around the common intersection.
#
plot(c(-3,3) + lambda, c(-3,3) + alpha, type="n", asp=1,
     main=paste(signif(beta,3), collapse=", "),
     xlab=expression(x[1]), ylab=expression(y))
#
# Plot the model for various values of x[2].
#
invisible(sapply(seq(-2,2,length.out=8)/beta[4], function(x2) {
  curve(beta[1] + beta[2]*x + beta[3]*x2 + beta[4]*x*x2, add=TRUE)
}))
#
# Highlight the common point.
#
points(lambda, alpha, pch=21, bg="Red")

