We propose median polish for functional multivariate analysis of variance (FMANOVA) with the implementation of depth for multivariate functional data. As an alternative to classical mean estimation, functional median polish estimates the functional grand effect and factor effects based on functional medians in one-way and two-way additive FMANOVA models. Median polish estimates in FMANOVA are visually unbiased, independently of the choice of multivariate functional depth. The corresponding mean-based and rank-based tests are generalized to evaluate whether the functional medians in various levels of the factors are the same. Simulation studies illustrate the robustness of our functional median polish in various scenarios, compared with the results from classical FMANOVA fitted by means. The results are evaluated both marginally and jointly. Three environmental datasets are considered to illustrate that our median polish is robust against outliers in practical implementations. Functional boxplots and heatmaps are two ways of visualizing the functional factors, depending on whether the functional data are curves or images, respectively.