def METHOD_NAME(domain, f, fh, num_pts=100):
    x  = np.linspace(*domain, num=num_pts)
    y  = f(x)
    yh = np.array([fh(xi) for xi in x])
    max_norm_err = np.max(abs(yh - y))
    print("Maximum error on evaluation grid: {}".format(max_norm_err))
    fig, ax = plt.subplots(1, 1)
    ax.plot(x, y, label='f(x)')
    ax.plot(x, yh, '.', label='f1_h(x)')
    ax.grid(True)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    fig.legend()
    fig.show()