Setting up boundary conditions to solve PDEs using method of lines

Question

Objective: To add boundary/initial conditions (BCs/ICs) to a system of ODEs

I have used the method of lines to convert a system of PDEs into a system of ODEs. The ODEs themselves involve a lot of variables so I will present simple versions here

For the purpose of this simulation all I am interested in is the structure of the code to set up the BCs/ICs, so I have removed any correlations etc and just put constants

For reference, the actual system is the flow of a gas through a packed bed. For terminology, I am calling the initial conditions the constant conditions along the bed initially (the bed temperature, the fluid temperature, the density of the fluid in the bed) and the boundary conditions would be the mass flow into the bed, the temperature of the flow into the bed.

I have two functions. df, which is used to set up the non IC/BC ODEs (generally happy with the structure of this one):

    def df(t,x,*constants):
##the time derivative of x
a,b,c,d = constants  ## unpack constants

##setting up solution array
x = x.reshape(2,number_of_nodes)
dxdt = np.zeros_like(x)

## internal nodes
for j in range(2,2*n,2): ##start at index 2 (to avoid ICs/BCs), go to 2*n and increment by 2
    ##equation 1
    dxdt[0][j] = a*(x[1][j-1] - b*x[1][j])
    ## equation 2 
    dxdt[1][j] = c*(x[1][j-1] - d*x[1][j])
    
return dxdt.reshape(-1)

And initialize_system, which is supposed to set up the BCs/ICs (which is currently incorrect):

def initialize_system(*constants):
#sets the appropriate initial/boundary conditions

a,b,c,d = constants #unpack constants
x0 = np.zeros((2, number_of_nodes))

y01 = 1
y02 = 1

#set up initial boundary values
# eq 1
x0[0][0] =  a*(y01 - b*x0[1][0]) 
# eq 2
x0[1][0] = c * (y02 - d*x0[1][0]) 

for j in range(2, 2*n, 2):
    # eq 1
    x0[0][j] = a*(y0CO2 - b*x0[1][j])
    # eq 2
    x0[1][j] = b*(y0H2O - d*x0[1][j])

return x0.reshape(-1)

My question is:

How can I correctly set up the ICs/BCs in initialize_system here?

Edit: There are very limited examples of implementing the method of lines in python on the internet. I would also highly appreciate any resources on this