Individual participant data (IPD) meta-analysis projects obtain, harmonise, and synthesise original data from multiple studies. Many IPD meta-analyses of randomised trials are initiated to identify treatment effect modifiers at the individual level, thus requiring statistical modelling of interactions between treatment effect and participant-level covariates. Using a two-stage approach, the interaction is estimated in each trial separately and combined in a meta-analysis. In practice, two complications often arise with continuous outcomes: examining non-linear relationships for continuous covariates and dealing with multiple time-points. We propose a two-stage multivariate IPD meta-analysis approach that summarises non-linear treatment-covariate interaction functions at multiple time-points for continuous outcomes. A set-up phase is required to identify a small set of time-points; relevant knot positions for a spline function, at identical locations in each trial; and a common reference group for each covariate. Crucially, the multivariate approach can include participants or trials with missing outcomes at some time-points. In the first stage, restricted cubic spline functions are fitted and their interaction with each discrete time-point is estimated in each trial separately. In the second stage, the parameter estimates defining these multiple interaction functions are jointly synthesised in a multivariate random-effects meta-analysis model accounting for within-trial and across-trial correlation. These meta-analysis estimates define the summary non-linear interactions at each time-point, which can be displayed graphically alongside confidence intervals. The approach is illustrated using an IPD meta-analysis examining effect modifiers for exercise interventions in osteoarthritis, which shows evidence of non-linear relationships and small gains in precision by analysing all time-points jointly.