How to calculate uncertainty of y values generated by fit functions

Question

I am using a fit function to calculate values used by an application in a manner similar to below:

    import numpy as np
    from numpy import random

    x = range(10)
    y = random.standard_normal(10) 
    w = random.standard_normal(10)/10 
    w = 1/w

    p,cov = np.polynomial.polynomial.polyfit(x=x,y=y,deg=1,w=w,full=True) 
    fun = np.polynomial.polynomial.Polynomial(p)

    new_x = 20 
    new_y = fun(new_x) 
    #y_1_sigma_uncertainty = ???

Is there a way to use the covariance matrix to calculate an uncertainty associated with values calculated by fun? Is there another way to go about this? I have done quite a bit of searching, but I am probably not asking the question correctly. I am not a stats person so I am hoping my example is useful in clarifying what I am trying to ask.

Thanks, gl