In addition to the good answer from Todd B, here's a process you can follow, and some more of the background math. 

### A reasonable process to follow

 1. Fit a Cox model (`fit.unstrat`) with all of the covariates and no stratifying variables.

 2. Use `cox.zph(fit.unstrat)` to check for violations of the proportional hazards assumption. If you see e.g. `rx` has p<0.05, then `rx` violates the proportional hazards assumption and should not be included as a covariate in the Cox model. A reasonable follow-up is to try changing it to a stratifying variable, e.g., `fit.strat <- update(fit.unstrat, . ~ . - rx + strata(rx))`.

 3. Next, it's appropriate to evaluate whether multiple models improve over the stratified model. The stratified Cox model allows the two `rx` groups to have different baseline hazards / baseline survivals, but it still only calculates a single effect $\beta_\text{age}$ of age on survival time. You can think of this as being the mean effect of age averaged across both `rx` groups. This raises a question: do you get a statistically significantly better fit when you allow the two `rx` groups to have different covariate effects ($\beta_\text{age, rx=2}$ and $\beta_\text{age, rx=1}$)? To answer this question, you:
    - fit the stratified model, 
    - fit two separate Cox models for datasets subsetted by `rx`, and 
    - compare them in a likelihood ratio test, where the p-value is $P\left(\chi^2_{(k-1)(q)}\right)>\ell_\text{stratified model} - \sum \ell_\text{unstratified models}$, with 
$k$ being the number of stratifying groups (just two here--`rx`=1 and `rx`=2), $q$ being the number of covariates (just 1 here--$z_\text{age}$), and $\ell$ is a model's log-likelihood.

Here's some sample code for the likelihood ratio tests. In this example (pretending `rx` violates the proportional hazards assumption), we wouldn't choose the multiple models, because while they improve the log likelihood, it's not a statistically significant improvement in the fit (p=0.1575):

    library(survival)
    fit.strat <- coxph(Surv(futime, fustat) ~ age + strata(rx), data=ovarian)
    fit.grp1 <- coxph(Surv(futime, fustat) ~ age, subset=(rx==1), data=ovarian)
    fit.grp2 <- coxph(Surv(futime, fustat) ~ age, subset=(rx==2), data=ovarian)
    
    LL.strat <- fit.strat$loglik[2]
    LL.unstrat <- fit.grp1$loglik[2] + fit.grp2$loglik[2]
        
    X2 <- -2*(LL.strat - LL.unstrat)
    n_groups <- length( unique(ovarian$rx) )
    n_params <- length( fit.strat$coef )
    p_val <- 1 - pchisq(X2, df=(n_groups-1)*n_params)

### Background math
**Standard** Cox proportional hazards model:

$$h(t|Z) = h_0 \exp(\beta_\text{age}Z_\text{age} + \beta_\text{rx}Z_\text{rx})$$

    fit.unstrat <- coxph(Surv(futime, fustat) ~ age + rx, data=ovarian)

 - The quantity $h_0$ is called the baseline hazard. It corresponds with the baseline survival, which is the survival of the reference group whose covariates are all equal to 0.

 - For a Cox model, the baseline hazard is estimated from the data non-parametrically (without assuming any distribution).

 - The equation above implies:
     1. (hazard for `rx==2` group)/(hazard for `rx==1` group) is a constant. In other words, group 2's hazard is always higher/lower than group 1's hazard by a constant factor.

     2. The two `rx` groups have the same baseline hazard and baseline survival.

**Stratified** Cox proportional hazards model (you're actually fitting two equations that share a single $\beta_\text{age}$):

$$h_\text{rx=2}(t|Z) = h_{0,\text{rx=2}} \exp(\beta_\text{age}Z_\text{age}) \\ 
h_\text{rx=1}(t|Z) = h_{0,\text{rx=1}} \exp(\beta_\text{age}Z_\text{age})$$

    fit.strat <- coxph(Surv(futime, fustat) ~ age + strata(rx), data=ovarian)

These equations imply:
  1. The two rx groups have different baseline hazards / baseline survivals.
  2. The two rx groups share a *single factor*, $\exp(\beta_\text{age})$, by which the hazard is always higher/lower for a 1-unit increase in age. In other words, we're specifying that both `rx` groups have the same $\beta_\text{age}$.

**Multiple models** (this is just a term I'm making up):

$$h_\text{rx=2}(t|Z) = h_{0,\text{rx=2}} \exp(\beta_\text{age, rx=2}Z_\text{age}) \\ 
h_\text{rx=1}(t|Z) = h_{0,\text{rx=1}} \exp(\beta_\text{age, rx=1}Z_\text{age})$$

    fit.grp1 <- coxph(Surv(futime, fustat) ~ age, subset=(rx==1), data=ovarian)
    fit.grp2 <- coxph(Surv(futime, fustat) ~ age, subset=(rx==2), data=ovarian)

### Short note about time-dependent covariates

In addition to fitting a stratified model, another typical way to remediate a violation of the proportional hazards assumption is to (1) convert the violating variable into a time-dependent variable, and then (2) fit a new Cox PH model and check if the time-dependent covariate satisfies the PH assumption. Time-dependent covariates are an entirely different topic, but they can get invoked for the same reason as a stratified Cox model. Picking between these two often depends on context. If you think the two rx groups truly have different baseline survivals, then a stratified model makes sense. If there's a time component to `rx` (e.g., the data also contain a variable indicating how much time passed before they started treatment, or we think `rx` might have a cumulative effect such that the treatment only provides a benefit after some amount of time has passed), then a time-dependent variable might make more sense.