Q:

Is there a way to align 2 independent axes in an R graph?

I'm plotting a graph with different axis. The problem is that I want the 2 axes to be crossing one point, the rest doesn't really matter. Is it possible?
Here is a reproducible code: 
plot(x = -10:10, y = -10:10)
abline(v=0,lty = 2)
par(new =TRUE)
plot(x = -10:50, y = seq(-5,5,length.out = length(-10:50)), xaxt = "n", yaxt = "n", bty ="n")
abline(v=0,lty = 3)
axis(3, col="red",col.axis="red",las=2, cex.axis = 1)
axis(4, col="red",col.axis="red",las=2, cex.axis = 1)

Here is the output: 

I basically want the vertical lines to be crossing each other at 0. 
Is there another way to write this (which is not really accurate, I want something that can be done automatically, other than setting a xlim manually, since I have a lot of them to do. Also, when the image is resized, it's losing the alignment with the 2 0s): 
plot(x = -10:10, y = -10:10)
abline(v=0,lty = 2)
par(new =TRUE)
plot(x = -10:50, y = seq(-5,5,length.out = length(-10:50)), 
     xaxt = "n", yaxt = "n", bty ="n",
     xlim = c(-50,50))
abline(v=0,lty = 3, lwd = 5)
axis(3, col="red",col.axis="red",las=2, cex.axis = 1)
axis(4, col="red",col.axis="red",las=2, cex.axis = 1)

The output should use something similar to what biplot.prcomp is using to align the arrows and the axis: 

EDITED
With a PCA, it's only working for the y axis, not the x axis.
new_lim <- function(a, type = 1) {
  newdata_ratio <-  NULL
  i <- type * 2 - 1
  old_lim <- par("usr")[i:(i+1)] + c(diff(par("usr")[i:(i+1)]) * 0.04 / 1.08, 
                                     diff(par("usr")[i:(i+1)]) * -0.04 / 1.08)
  old_ratio <- old_lim[1] / old_lim[2]
  newdata_ratio <- if (max(a) <= 0) -1.0e+6 else min(a) / max(a)
  if (old_ratio >= newdata_ratio ) {
    new_min <- min(a)
    new_max <- min(a) / old_ratio
  } else {
    new_min <- max(a) * old_ratio
    new_max <- max(a)
  }
  c(new_min, new_max)
}
s1= rnorm(50,mean = 12)
s2= rnorm(50, mean = 17)
s3= rnorm(50, mean = 20)
library(vegan)
pca=rda(cbind(s1,s2,s3))
pca.scoop=scores(pca, scaling = 2)
biplot(pca)

par(mar=c(4, 4, 4, 4))
plot(pca, xlab = "x1", ylab = "y1",
     type = c("p"),
     main= "main",
     scaling = 2,
     choices = c(1,2),
     xlim =c(min(pca.scoop$sites[,1]),max(pca.scoop$sites[,1])),
     ylim = c(min(pca.scoop$sites[,2]),max(pca.scoop$sites[,2])),
     bty = "o",#"l"
     pch=4)
abline(v = 0, lty = 2); abline(h = 0, lty = 2)

x2 <- -10:20
y2 <- seq(40, 10, length.out = length(x2))
par(new =TRUE)
plot(x2, y2, 
     xlim = new_lim(x2), 
     ylim = new_lim(y2, 2), axes = F, ann = F)
axis(3, col = "red", col.axis = "red") # axes=F is equivalent to xaxt="n", yaxt="n" and bty="n"
axis(4, col = "red", col.axis = "red") # ann=F is equivalent to xlab=NA and ylab=NA
mtext("x2", side = 3, line = 2.5, col = "red")
mtext("y2", side = 4, line = 2.5, col = "red")
# box(bty="7", col="red")  # if you want.

A:

Is it possible to align the second graph based on the first one?

Yes, it is. But it's a litte complex because of needing logical judgment. 
This function, new_lim(a, type), calculates first graph's xlim (or ylim) and ratio of minus to plus from graphic parameters. And it judges which value it uses as second graph's parameter, min(second data) or max, and calculate another value from the first graph's ratio. The output are second graph's min and max value, _lim. The argument a is a second x or y data. type = 1 (default; omittable) is for xlim, type = 2 (type = is omittable) is for ylim.
new_lim <- function(a, type = 1) {
  newdata_ratio <-  NULL
  i <- type * 2 - 1
  old_lim <- par("usr")[i:(i+1)] + c(diff(par("usr")[i:(i+1)]) * 0.04 / 1.08, 
                                     diff(par("usr")[i:(i+1)]) * -0.04 / 1.08)
  old_ratio <- old_lim[1] / old_lim[2]
  newdata_ratio <- if (max(a) <= 0) -1.0e+6 else min(a) / max(a)
  if (old_ratio >= newdata_ratio ) {
    new_min <- min(a)
    new_max <- min(a) / old_ratio
  } else {
    new_min <- max(a) * old_ratio
    new_max <- max(a)
  }
  c(new_min, new_max)
}

[Note]
This function needs to exist of first graph and include zero of first data's range. It is no problem that second data's range doesn't include zero.
x2 <- -40:20
y2 <- seq(40, 10, length.out = length(-40:20))

par(mar=c(4, 4, 4, 4))
plot(x = -15:5, y = -5:15, xlab = "x1", ylab = "y1")
abline(v = 0, lty = 2); abline(h = 0, lty = 2)

par(new =TRUE)
plot(x2, y2, xlim = new_lim(x2), ylim = new_lim(y2, 2), axes = F, ann = F)
axis(3, col = "red", col.axis = "red") # axes=F is equivalent to xaxt="n", yaxt="n" and bty="n"
axis(4, col = "red", col.axis = "red") # ann=F is equivalent to xlab=NA and ylab=NA
mtext("x2", side = 3, line = 2.5, col = "red")
mtext("y2", side = 4, line = 2.5, col = "red")
# box(bty="7", col="red")  # if you want.

EDITED  [CAUTION]
When you use this function with plot(class.rda) and change the aspect by Rsutdio, the output becomes different from what you want.
x2 <- -10:20
y2 <- seq(40, 10, length.out = length(x2))
library(vegan) 
s1= rnorm(50,mean = 12); s2= rnorm(50, mean = 17); s3= rnorm(50, mean = 20)
pca=rda(cbind(s1,s2,s3))
pca.scoop=scores(pca, scaling = 2)
biplot(pca)
par(new =TRUE)
plot(x2, y2, 
     xlim = new_lim(x2), 
     ylim = new_lim(y2, 2), axes = F, ann = F)
axis(3, col = "red", col.axis = "red")
axis(4, col = "red", col.axis = "red")
mtext("x2", side = 3, line = 2, col = "red")
mtext("y2", side = 4, line = 2, col = "red")

