result = service.job(<some zne job id>).result()
pub_result = result[0]

nfs = result.metadata["resilience"]["zne"]["noise_factors"]
nfs_x = result.metadata["resilience"]["zne"]["extrapolated_noise_factors"]

models = result.metadata["resilience"]["zne"]["extrapolator"]
plt.figure(figsize=(5 * len(models), 4))

for idx_model, model in enumerate(models):
    plt.subplot(1, len(models), idx_model + 1)
    for idx in range(pub_result.data.evs_noise_factors.shape[2]):
        lines = plt.errorbar(
            nfs,
            pub_result.data.evs_noise_factors[0,0,idx],
            fmt='.',
            yerr=pub_result.data.stds_noise_factors[0,0,idx]
        )

        plt.plot(
            nfs_x,
            means:=pub_result.data.evs_extrapolated[0,0,idx,idx_model,:],
            color=(color:=lines[0].get_color()),
            alpha=0.3
        )

        plt.fill_between(
            nfs_x,
            means-(stds:=pub_result.data.stds_extrapolated[0,0,idx, idx_model, :]),
            means+ stds,
            color=color,
            alpha=0.1
        )

        plt.ylim([-0.05,1.05])