My teacher gave the R language as follows
p1=130
p2=113
set1=seq(2,402,10)
result=matrix(0,1,8)
for (n1 in set1){
rv1=matrix(0,100,1)
rv2=matrix(0,100,1)
for (j in 1:100)  {

x=mvrnorm(n=n1,rep(0,p1),diag(1,p1))
y=mvrnorm(n=n1,rep(0,p2),diag(1,p2))

num1=sum((t(y)%*%x)^2)
den1=sqrt(sum((x%*%t(x))^2)%*%sum((y%*%t(y))^2))
V1=(x%*%t(x))-diag(diag(x%*%t(x)))
W1=(y%*%t(y))-diag(diag(y%*%t(y)))
V2=matrix(V1,length(V1),1)
W2=matrix(W1,length(W1),1)
num2=t(V2)%*%W2
den2=sqrt((t(V2)%*%V2)%*%(t(W2)%*%W2))

rv1[j]=num1/den1
rv2[j]=num2/den2
}
rv3=(p1%*%p2)/sqrt((p1^2+(n1+1)%*%p1)*(p2^2+(n1+1)%*%p2))
result=rbind(result,cbind(n1,mean(rv1),min(rv1),max(rv1),rv3,mean(rv2),min(rv2),max(rv2)))
}
result=result[-1,]
result
cbind(seq(2,402,10),result[,5]-result[,2])
x1=result[,1]

plot(x1, result[,2], type='l', col="red")
points(x1, result[,3], type='l', col="red")
points(x1, result[,4], type='l', col="red")
plot(x1, result[,6], type='l',ylim=c(-1,1), xlim=c(0,500), col="red")
points(x1, result[,7], type='l',ylim=c(-1,1), xlim=c(0,500), col="red")
points(x1, result[,8], type='l',ylim=c(-1,1), xlim=c(0,500), col="red")