human
I am working on the file "LDL_transport_model.ipynb".
The first task is:
Assembling and solving the discrete system
The LDL_Sim class is responsible for solve the differential equation. 
In the first step the simulation region is discretized by the method discretize. In that function flux continuity between layers and boundary conditions are implemented. It is done in a following way: 

the layers are encoded as list layers in discretize, containing indices of boundaries for given discretization

system parameters depending on space ($k,\vec u,\sigma,D_{eff}$) are sampled

the diagonals of matrix are assembled using finite differences in space, neglecting for a moment the space dependence 

at regions boudaries the equation is replaced with 
  $$J_L = J_R$$
  it takes a form for a boundary at index j:
  $$ c_j v_{j-1} - D_{j-1}\frac{c_j-c_{j-1}}{dx} =  c_j v_{j+1} - D_{j+1}\frac{c_{j+1}-c_{j}}{dx}  $$
  collecting terms with the same $c$:
  $$ c_j \left( v_{j-1}- v_{j+1}  - \frac{D_{j-1}+D_{j+1}}{dx}\right) + c_{j-1}\frac{D_{j-1}}{dx} + c_{j+1}\frac{D_{j+1}}{dx} = 0
    $$

therefore we need to enforce the above equation on the system for each boundary inside the domain

the system matrix is stored as scipy.sparse dia_matrix

The plot_c function is used for plot the concentration profiles.
Can you write Python code for it?
 Can you provide a the Python code you requested?
assistant
#Here is the Python code you requested.

class LDL_Sim(object):
   
    def __init__(self, pars):
        self.pars = pars
        self.c_st = None
    def discretize(self,N=2000):
        self.N = N
        k = np.ones(N)
        v = np.ones(N)
        Dyf = np.ones(N)

        D,V,sigma,L,k_react = self.pars.get_params()
        l = np.sum(L) 
        self.l = l
        self.x=np.linspace(0,l,N)
        
        layers=[0]+list( np.ceil( (N*(np.cumsum(L)/sum(L)))).astype(np.int32) )
    
        for i,(l1,l2) in enumerate(zip(layers[:],layers[1:])):
            k[l1:l2] = k_react[i]
            v[l1:l2] = (1.0-sigma[i])*V[i]
            Dyf[l1:l2] = D[i]
        dx2_1 = (N-1)**2/l**2
        dx_1 = (N-1)/l

        diag_l = np.ones(N)*(np.roll(Dyf,-1)*dx2_1)
        diag   = np.ones(N)*(-2.*Dyf*dx2_1 - k + v*dx_1)
        diag_u = np.ones(N)*(np.roll(Dyf,1)*dx2_1 - np.roll(v,1)*dx_1)

        # Layer's junctions
        for j in layers[1:-1]:
            diag[j] = v[j-1]-v[j+1]-(Dyf[j-1]+Dyf[j+1])*dx_1
            diag_l[j-1] = Dyf[j-1]*dx_1
            diag_u[j+1] = Dyf[j+1]*dx_1
        #Boundary Conditions
        diag[0] = 1
        diag[-1] = 1
        diag_u[0+1] = 0
        diag_l[0-2] = 0
   
        self.L = dia_matrix((np.array([diag_l,diag,diag_u]),np.array([-1,0,1])), shape=(N,N))
        
    def solve_stationary(self,bc=[1,0]):
        b = np.zeros(self.N)
        b[0],b[-1] = bc
        L = self.L.tocsr()
        self.c_st = sp.sparse.linalg.linsolve.spsolve(L,b)

        
    def plot_c(self,yrange=(0,0.2),xrange=(0,214),filename=None, color='red', alpha=0.2, style='-'):
        
        i1,i2 = int(xrange[0]/self.l*self.N),int(xrange[1]/self.l*self.N)
        plt.plot(self.x[i1:i2],self.c_st[i1:i2],color=color,linewidth=2, ls=style)
        plt.ylim( *yrange)
        plt.xlim( *xrange)
        L=self.pars.L
        d=[0]+np.cumsum(self.pars.L).tolist()
        colors=['m','g','b','w']
        for i,(l1,l2) in enumerate(zip(d[:],d[1:])):
            
            plt.bar([l1,],yrange[1],l2-l1, color=colors[i], linewidth=0.3, ...