dist = "weibullPH"
rather than dist = "weibull"
provides scale and shape parameters that can be used as λ and γ, respectively in the parameterisation below to calculate the survival function. (k AKA γ).
Using larynx data:
library("KMsurv")
library("survival")
library("flexsurvreg")
data(larynx)
larynx$late_stage <- factor(larynx$stage == 3 | larynx$stage == 4) #create a binary variable for illustration
swei <- flexsurvreg(s ~ late_stage,dist='weibullPH', data=larynx)
plot(swei, type = "survival", xlim = c(0,10)) #plot the output from the WeibullPH model
#Now extract the lambda and gamma from swei$res
lambda <- swei$res[2]
gamma <- swei$res[1]
t<- larynx$time
You can then plot the calculated survival function for group 1 using $ S(t) = exp(-\lambda t^{\gamma}) $ and for group 2 by introducing the effect of the coefficient through raising the previous equation to the power of the linear predictor $ S(t) = exp(-\lambda t^{\gamma})^{exp(\hat{\beta}x_{i}) } $
group1 <- exp(-lambda*(t^gamma))
group2 <- group1^(exp(swei$coefficients[3]))
plot(swei, type = "survival", xlim = c(0,10))
par(new=TRUE)
plot(t, group1 , ylim=c(0,1), xlim = c(0,10))
par(new=TRUE)
plot(t, group2, ylim=c(0,1), xlim = c(0,10))

You can see that your manual survival plot overlies that of the function perfectly.
You could achieve the same thing using using the scale and shape values provided with dist="Weibull" using the alternative parameterisation $ S(t) = exp(-(\frac{t}{\lambda})^{\gamma }) $