par(mfrow=c(1,1),cex=1.3)
n=100
R0=c(-0.99,-.95,-0.8,-0.6,-0.2,0.1,0.3,0.5,0.7,0.9)
m=sample(1:10,1)
R=R0[m]

# generate a sample whose sample correlation is guaranteed equal to R
x=rnorm(n)
y=rnorm(n)
sxx=var(x)
syy=var(y)
sxy=sqrt(sxx*syy)*cor(x,y)
f1=function(a){a*sxx+sqrt(1-a^2)*sxy-R*sqrt(a^2*sxx+(1-a^2)*syy+2*a*sqrt(1-a^2)*sxy)*sqrt(sxx)}
library(rootSolve)
a0=uniroot(f1,c(-1,1))$root
z=a0*x+sqrt(1-a0^2)*y
plot(x,z,pch=20)