How to plot PDF where a parameter is a function of another variable

They are plotting the 5, 50, and 95 percentiles of the "cumulative proportion failing function", i.e. $\text{Pr}(Y_i\le y)=F(y;\mu(x_i),\sigma(x_i))$, for which you've provided the definition above.

For a fixed $x$ that becomes a relatively simple cumulative normal density function with a single unknown, $y$, that will determine the actual proportion. Solving for the desired percentile will give you the lines in the plot.

I'll use the parameter estimates from the paper to calculate $\mu(x_i)$ and $\sigma(x_i)$ since I understood you were able to reproduce these already. I'm more of an R person but this should be very straightforward in any statistical software:

## Parameter estimates from the paper
B0m <- 14.75
B1m <- -1.39
B0s <- 10.97
B1s <- -2.50
G <- 75.71

## mu(x) and sigma(x)
mufn <- \(x) ifelse(x > G, B0m + B1m*log(x-G), NA)
sigfn <- \(x) exp(B0s + B1s*log(x))

## cumulative proportion failing function F(y; mu(x), sigma(x))
Ffn <- \(y, x) pnorm((log(y) - mufn(x)) / sigfn(x))

## Wrapper to solve above for a given percentile
optfn <- \(y, x, t=0.05) Ffn(y, x) - t

## Fixed x-es to solve at
x <- seq(76, 145, length.out=1E2)
## Interval to search over for root
intr <- exp(c(-10, 30))

## 5% line
lower <- vapply(x, \(x) uniroot(optfn, intr, x=x, t=0.05)$root, numeric(1))
## 50% line
middle <- vapply(x, \(x) uniroot(optfn, intr, x=x, t=0.50)$root, numeric(1))
## 95% line
upper <- vapply(x, \(x) uniroot(optfn, intr, x=x, t=0.95)$root, numeric(1))

Plotting that looks convincingly similar to the paper's figure:

Raw data

fatigue <- tibble::tribble(
  ~stress, ~cycles, ~failed,
   80.3, 211629, 1,
   80.6, 200027, 1,
   80.8,  57923, 0,
   84.3, 155000, 1,
   85.2,  13949, 1,
   85.6, 112968, 0,
   85.8, 152680, 1,
   86.4, 156725, 1,
   86.7, 138114, 0,
   87.2,  56723, 1,
   87.3, 121075, 1,
   89.7, 122372, 0,
   91.3, 112002, 1,
   99.8,  43331, 1,
  100.1,  12076, 1,
  100.5,  13181, 1,
  113.0,  18067, 1,
  114.8,  21300, 1,
  116.4,  15616, 1,
  118.0,  13030, 1,
  118.4,   8489, 1,
  118.6,  12434, 1,
  120.4,   9750, 1,
  142.5,  11865, 1,
  144.5,   6705, 1,
  145.9,   5733, 1
)

Plot code

plot(NULL, NULL, xlim=c(5, 500), ylim=c(78, 150),
     log="x", xlab="kilocycles", ylab="stress")
points(fatigue$cycles/1000, fatigue$stress,
       pch=ifelse(fatigue$failed, 19, 17),
       col=ifelse(fatigue$failed, "black", "red"), cex=0.8)
lines(lower/1000, x, lty="36")
lines(middle/1000, x)
lines(upper/1000, x, lty="36")