code: def coregionalization_toy(optimize=True, plot=True):
      
      #build a design matrix with a column of integers indicating the output
      X1 = np.random.rand(50, 1) * 8
      X2 = np.random.rand(30, 1) * 5
  
      #build a suitable set of observed variables
      Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05
      Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2.
  
      m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2])
  
      if optimize:
          m.optimize('bfgs', max_iters=100)
  
      if plot:
          slices = GPy.util.multioutput.get_slices([X1,X2])
          m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0})
          m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca())
      return m
docs: A simple demonstration of coregionalization on two sinusoidal functions.