We extend cokriging analysis and multivariable spatial prediction to the case where the observations at each sampling location consist of samples of random functions, that is, we extend two classical multivariable geostatistical methods to the functional context. Our cokriging method predicts one variable at a time as in a classical multivariable sense, but considering as auxiliary information curves instead of vectors. We also propose an extension of multivariable kriging to the functional context by defining a predictor of a whole curve based on samples of curves located at a neighborhood of the prediction site. In both cases a non-parametric approach based on basis function expansion is used to estimate the parameters, and we prove that both proposals coincide when using such an approach. A linear model of coregionalization is used to define the spatial dependence among the coeficients of the basis functions, and therefore for estimating the functional parameters. As an illustration the methodological proposals are applied to analyze two real data sets corresponding to average daily temperatures measured at 35 weather stations located in the Canadian Maritime Provinces, and penetration resistance data collected at 32 sampling sites of an experimental plot.