You can do this "by hand". Here is one possibility using the example data from coxph(). Essentially, you need to make a dataset that holds constant all model variables at values you want to use and then vary time from the minimum to maximum in the data (here 0:4). Then, you can predict the survival probabilities for these data and plot.
library(ggplot2)
library(survival)
test1 <- list(time=c(4,3,1,1,2,2,3),
status=c(1,1,1,0,1,1,0),
x=c(0,2,1,1,1,0,0),
sex=c(0,0,0,0,1,1,1))
# Fit a stratified model
m1 <- coxph(Surv(time, status) ~ x + strata(sex), test1)
tmp1 <- data.frame(
time = 0:4,
status=0,
x = 1,
sex = 0)
tmp1$fit <- predict(m1, newdata=tmp1, type="survival")
ggplot(tmp1, aes(x=time, y=fit)) +
geom_line() +
geom_point() +
theme_classic()

Below, we pick a different value of x, generate predictions and put the data together with the prediction data above. This could be done in a single step if you wanted.
tmp2 <- data.frame(
time = 0:4,
status=0,
x = 2,
sex = 0)
tmp2$fit <- predict(m1, newdata=tmp2, type="survival")
tmp <- rbind(tmp1, tmp2)
Now, we could plot both sets of predictions on the same plot.
ggplot(tmp, aes(x=time, y=fit, colour=as.factor(x))) +
geom_line() +
geom_point() +
theme_classic() +
labs(colour="X")

Created on 2022-04-29 by the reprex package (v2.0.1)