Simulation of random variables using Polar method

I have the following algorithm

Step 1. Generate u1 and u2~U(0,1)

Step 2. Define v1=2u1-1, v2=2u2-1 and s=v1^2+v2^2

Step 3. If s>1, come back to Step 1.

Step 4. If s<=1, x=v1(-2logs/s)^(1/2) and y=v2(-2logs/s)^(1/2)

Here is my approach to implement this algorithm in R:

    PolarMethod1<-function(N)
{

  x<-numeric(N)
  y<-numeric(N)
  z<-numeric(N)

  i<-1

  while(i<=N)
  {u1<-runif(1)
  u2<-runif(1)
  v1<-(2*u1)-1
  v2<-(2*u2)-1
  s<-(v1^2)+(v2^2)

  if(s<=1)
  {
    x[i]<-((-2*log(s)/s)^(1/2))*v1
    y[i]<-((-2*log(s)/s)^(1/2))*v2
    z[i]<-(x[i]+y[i])/sqrt(2) #standarization
    i<-i+1
  }
  else
    i<-i-1
  }

  return(z)
}
z<-PolarMethod1(100000)
hist(z,freq=F,nclass=10,ylab="Density",col="purple",xlab=" z values")
curve(dnorm(x),from=-3,to=3,add=TRUE)

The code, fortunately, does not mark any error and works quite well with N=1000 but when I change to N=10000, instead of making a better approach to the curve displays:

contrast with N=1000 displays:

Why is that?

Is there something wrong with my code? It's supposed to be better adjusted when N increases.

Note:I added the z in the code to include both variables in the output.