There might be some small peformance improvements you can make (ex. I don’t think you need the for loop in your computation of mu), but in general GPs become unfeasible computationally with more than 100 or so sample points. You might want to check out splines as an approximation for such cases.