Matplotlib animated subplots using multiple lines

I'm trying to compare several numerical schemes with exact solution of a PDE. To do that, I'm trying to make an animation, in which there are multiple subplots and each subplot has two functions evolving. (numerical and exact)

I found out how to make multiple lines animating here, Matplotlib multiple animate multiple lines

And I also found out how to make multiple subplots. animated subplots using matplotlib

But I'm not sure how to do both. They looks so confusing when looked together (I'm not that familiar with matplotlib, especially FuncAnimation) And also I want to make a subplot number as variable, so user can input any number of functions.

my code is this. Currently this code generates plot for each numerical scheme, but I want it to make only one plot. So the 'plot' function would be called only once,

plot(*arguments, FTCS, upwind, anotherfunc, anotherfunc2, ...)

this way.

Sorry for the messy code. I attached all of them just to make a minimum working example. Only 'plot' function is relevant to my question.

*edited code.

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



#dt=1/5
dx=1/400
tm=9
xm=1
C=0.8
u=.1
dt=dx*C/u
arguments=dt,dx,tm,xm,C

def plot(dt,dx,tm,xm,C,*funcs):
    plt.figure()
    fig,axes=plt.subplots(len(funcs))
    i=1
    for ax,func in zip(axes, funcs):
    #    plt.subplot(len(funcs),1,i)
        x1,T1=func(dt,dx,tm,xm,C)
        plotlays, plotcols = [2], ["black","red"]
        lines = []
        ax=plt.axes(xlim=(0,1),ylim=(-1.1,1.1))
        line,=ax.plot([],[],lw=2)
        for index in range(2):
            lobj = ax.plot([],[],lw=2,color=plotcols[index])[0]
            lines.append(lobj)

        def init():
            for line in lines:
                line.set_data([],[])
            return lines

        def animate(i):
            x=x1
            y=T1[i]
            y2=T2[i]
            xlist=[x,x]
            ylist=[y,y2]

            for lnum,line in enumerate(lines):
                line.set_data(xlist[lnum],ylist[lnum])
            return lines

        anim=animation.FuncAnimation(fig,animate,init_func=init,frames=np.shape(T1)[0],interval=20,blit=False)

        plt.title('%s' % func.__name__)
        plt.xlabel('x')
        plt.ylabel('T')
        i+=1

    plt.show()



    return



def finit(x):
    if 0<=x<=.1:
        return np.sin(10*np.pi*x)
    else:
        return 0

def FTCS(dt,dx,tm,xm,C):
    x=np.arange(0,xm+dx,dx)
    t=np.arange(0,tm+dt,dt)
    lx=len(x)
    lt=len(t)
    T=np.zeros([lt,lx])
    for i in range(0,lx):
        T[0,i]=finit(x[i])
    for n in range(0,lt-1):
        for j in range(1,lx-1):
            T[n+1,j]=T[n,j]-.5*C*(T[n,j+1]-T[n,j-1])
    return x,T

def upwind(dt,dx,tm,xm,C):
    x=np.arange(0,xm+dx,dx)
    t=np.arange(0,tm+dt,dt)
    lx=len(x)
    lt=len(t)
    T=np.zeros([lt,lx])
    for i in range(0,lx):
        T[0,i]=finit(x[i])
    for n in range(0,lt-1):
        for j in range(1,lx-1):
            T[n+1,j]=(1-C)*T[n,j]+C*T[n,j-1]
    return x,T

def exact(dt,dx,tm,xm,C):
    x=np.arange(0,xm+dx,dx)
    t=np.arange(0,tm+dt,dt)
    u=C*dx/dt
    T=np.zeros([len(t),len(x)])
    for n in range(0,len(t)):
        for j in range(0,len(x)):
            if u*t[n]<x[j]<=u*t[n]+.1:
                T[n,j]=np.sin(10*np.pi*(x[j]-u*t[n]))
            else: T[n,j]=0
    return T

plot(*arguments,FTCS,upwind)