The R competing risks vignette shows ways to calculate both Cox (Section 3) and Fine-Gray (Section 4) models in a competing-risks scenario.* The Fine-Gray model performs a Cox regression on a specially formatted and weighted data set that allows an individual to be incorporated into the risk set for an event type even after experiencing a different type of event.
As the two model types are both based on Cox proportional hazards regression, predictions from the two model types are done similarly, as also explained in the vignette. The hazard ratios aren't sufficient to compute cumulative incidence in either case, as they operate around a baseline hazard (resp. subdistribution hazard) that must be estimated from the data and those hazard ratios. For a Cox model with competing risks, it's important to start with a model that combines all competing risks so that the cumulative incidence can be estimated properly. In R, as shown in the vignette, the survfit() function accepts either a multi-state Cox model or a Fine-Gray model to generate the cumulative incidence over time for specified covariate values, which can be plotted.
That's one interpretation of what you call a "cumulative incidence function for each of the independent variables." If you want to plot incidence at some specific time as a function of covariate values, instead of as a function of time for specified covariate values, you will have to do more work; I don't think that the survival package directly supports predictions for multi-state Cox models at specific times with a predict() function. Near the end of Section 3, the vignette shows an example of how to extract probabilities of each state at a specific time from a summary() of a survfit() result, which could be used for that purpose.
*I wasn't able to get a copy of the reference to which you linked, so I can't say exactly how the analysis was done in that case.
