With the data, this will be difficult because the axes differ by 44 orders or magnitude
suppressPackageStartupMessages({
library(ggplot2)
})
# constants
L <- 1.210^(-8)
u0 <- 1.610^(-19) * 23210^(-3)
m_star <- 9.10910^(-31)
hbar <- (6.6310^(-34)) / (2 * pi)
me <- 0.0679 * 10910^(-31)
mh <- 0.59 * 10910^(-31)
theta0 <- sqrt(((m_star) * (u0) * (L^(2))) / (2 * (hbar^(2))))
generate_theta <- function(k) (k*L)/2
runs <- seq(1,100)
dat <- data.frame(x = generate_theta(runs),y = theta0*cos(generate_theta(runs)))
plot(dat$x,dat$y)

p <- ggplot(dat)
p + geom_line(aes(x,y)) +
theme_minimal()

p + geom_line(aes(runs,x)) +
theme_minimal()
