Great! It’s not working because the GP class takes the covariance function, not the covariance matrix, which is what you have after evaluating it at x. The reason is that sample_gp needs to re-evaluate the covariance function at the x values you want to predict at.
The best solution would be to write your own subclass of Covariance for this. See Linear for a good example. You would need to implement __init__ and full. If you have any questions or issues doing so, don’t hesitate to ask. In the meantime, I’ll see if we can add something to the library that would do the trick here.