How to combine multiple lines in a single legend entry for a quantity which depends on multiple variables?

Question

Goal

• Calculate two quantities q1 and q2, which both depend on two variables v1 and v2.
• Plot both quantities for a select number of values of variable 1, as function of variable 2.
• For each selected value for variable 1: list the value, as well as the average value for quantities 1 and 2 (+ error margins) in one legend entry

I have looked at this question, which shows desired output for a single (line + filled area), but could not get it to work with this more extensive problem.

Question

The main question is to create the legend as described above. Thanks!

Code

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

def calc_quant1(x,y):
    q1 = x + y
    q1err = np.sqrt(x+y)
    return q1,q1err

def calc_quant2(x,y):
    q2 = x - y
    q2err = np.sqrt((x+y)/2)
    return q2,q2err


# variables
var1 = np.linspace(-6,10,5)
var2 = np.linspace(10,20,21)

# resulting quantities (preallocation)
quant1 = np.zeros((len(var1),len(var2)))
quant2 = np.zeros((len(var1),len(var2)))
quant1err = np.zeros((len(var1),len(var2)))
quant2err = np.zeros((len(var1),len(var2)))
dy_dx = np.zeros(len(var2))

# plot colors
colors = cm.jet(np.linspace(0,1,len(var1)))

# reference value
xref = 12

# initialize figure
fig = plt.figure(figsize=(8,6))
fig.suptitle('Title')
ax = fig.add_subplot(111)
ax.axvline(x=xref,color='black',linestyle='dashed',label='$x={}$'.format(xref))
    
# compute and plot
for i in np.arange(len(var1)):
    
    # select "measured" voltage
    v1 = var1[i] * np.ones(len(var2))
    v1err = var1err[i] * np.ones(len(var2))
    
    # compute heat flux
    quant1[i,:], quant1err[i,:] = calc_quant1(v1,var2)
    quant2[i,:], quant2err[i,:] = calc_quant2(v1,var2)
    
    # compute gradient near reference point (suboptimal)
    x1 = var2[var2<xref][-1]
    x2 = var2[var2>xref][0]
    y1 = quant2[i][var2<xref][-1]
    y2 = quant2[i][var2>xref][0]
    dy_dx[i] = (y2-y1)/(x2-x1)
    
    # plot results
    label1 = r'$v_{{1}}={:.0f}$ [unit]'.format(var1[i])
    label2 = r'$q_{{1}}(10)={:.0f}$ [unit]'.format(quant1[i,0])
    ax.plot(var2,quant1[i,:], linestyle='dotted', linewidth=2, color=colors[i])
    ax.plot(var2,quant2[i,:], linestyle='solid', linewidth=1, color=colors[i])
    ax.fill_between(var2,quant1[i,:]-quant1err[i,:],quant1[i,:]+quant1err[i,:], alpha=0.16, color=colors[i])
    ax.fill_between(var2,quant2[i,:]-quant2err[i,:],quant2[i,:]+quant2err[i,:], alpha=0.16, color=colors[i],\
                    label=label1+'\n'+label2)


# finalize figure
ax.set_xlim(np.min(var2),np.max(var2))
ax.set_xlabel('$v_{2}$ [unit]')
ax.set_ylabel('$q_{2}$ [unit]')
plt.tight_layout()
ax.annotate(text=r'$Q_{1}=v_{1}+v_{2}$',
            xy=(0.50,0.77), xycoords='axes fraction',
            xytext=(1.05,0.95), textcoords='axes fraction',
            arrowprops=dict(arrowstyle='->',connectionstyle='arc3,rad=0.2'),
            fontsize=20)
ax.annotate(text=r'$Q_{2}=v_{1}-v_{2}$',
            xy=(0.70,0.35), xycoords='axes fraction',
            xytext=(1.05,0.85), textcoords='axes fraction',
            arrowprops=dict(arrowstyle='->',connectionstyle='arc3,rad=0.1'),
            fontsize=20)
leg = ax.legend(loc='upper right', bbox_to_anchor=(1.475,0.775), labelspacing=1)
ax.grid()
plt.show()

Output