I've just seen this question made by myself one year ago that I managed to deal with recently, I hope correctly. I put the script in case that it can be of help to another person and to see if somebody detects something wrong.
Following the example given in the question we can calculate the hazard ratio and their 95% interval with
age.spline = ns(flchain$age, df=3)
X.p = cbind(pdata$sex == 'M',
predict(age.spline, pdata$age),
predict(age.spline, pdata$age) * (pdata$sex == 'M'))
lHR = apply(matrix(rowSums( t(t(X.p) * coef(nfit3)) ), ncol = 2), 1, diff)
AGE.pred = (X.p[X.p[,1]==1,] - X.p[X.p[,1]==0,])
s2.age = apply(AGE.pred, 1, function(x){
t(x) %*% vcov(nfit3) %*% x
})
plot(50:99, exp(lHR), type = 'l', ylab = 'Hazard ratio', xlab = 'Age', lwd=2)
points(50:99, exp(lHR - 1.96 * sqrt(s2.age)), type= 'l', lty = 2, lwd=2)
points(50:99, exp(lHR + 1.96 * sqrt(s2.age)), type= 'l', lty = 2, lwd=2)
abline(h=1)
which ens up with the following plot

Extended after a question appeared in ResearchGate.
The approach can be further generalized to other packages. For example the package mgcv
, which allows the user to fit an additive cox model with certain smooth terms (p.e. thin plate regression splines).
library(mgcv)
nfittp = gam(futime~sex+s(age,by=sex,bs='tp'),
data=flchain, weights = death,
family = 'cox.ph', method = 'REML')
df.male = expand.grid(sex = 'M', age = 50:99)
df.female = expand.grid(sex = 'F', age = 50:99)
y0.male_ = predict(nfittp,newdata=df.male,type="lpmatrix")
y0.female_ = predict(nfittp,newdata=df.female,type="lpmatrix")
y0.male = predict(nfittp,newdata=df.male,type="link", se.fit=TRUE)
y0.female = predict(nfittp,newdata=df.female,type="link", se.fit=TRUE)
dplot = data.frame(
age = 50:99,
b.male = y0.male$fit - y0.female$fit,
s.male = sapply(1:length(50:99), function(i){
l = matrix(y0.male_[i,] - y0.female_[i,])
t(l) %*% vcov(nfittp) %*% l
})
)
par(mfrow=c(1,2))
plot(50:99, exp(y0.female$fit), type='l', log='y', lwd=2, xlab="Age", ylab="Relative risk", ylim=c(0.1,50))
points(50:99, exp(y0.female$fit + 1.96 * y0.female$se.fit), type='l', col=1, lwd=2, lty=2)
points(50:99, exp(y0.female$fit - 1.96 * y0.female$se.fit), type='l', col=1, lwd=2, lty=2)
points(50:99, exp(y0.male$fit), type='l', col='red', lwd=2)
points(50:99, exp(y0.male$fit + 1.96 * y0.male$se.fit), type='l', col=2, lwd=2, lty=2)
points(50:99, exp(y0.male$fit - 1.96 * y0.male$se.fit), type='l', col=2, lwd=2, lty=2)
legend(50, 40, c("Female", "Male"), lty=1, lwd=2, col=1:2, bty='n')
abline(h=1)
plot(50:99, exp(dplot$b.male), type = 'l', ylab = 'Hazard ratio', xlab = 'Age', lwd=2, ylim = c(0.8,2))
points(50:99, exp(dplot$b.male - 1.96 * sqrt(dplot$s.male)), type= 'l', lty = 2, lwd=2)
points(50:99, exp(dplot$b.male + 1.96 * sqrt(dplot$s.male)), type= 'l', lty = 2, lwd=2)
abline(h=1)
