This is my first post here and I hope that there is someone here who once had a similar problem.
I have successfully managed to visualise a DC membrane reactor using the solve_ivp solver. In other words, tube by tube and substance and temperature are exchanged via the interfaces. Now I was wondering what it would look like if I simulated the whole thing in countercurrent. My problems lie in the formulation of the BC and in the formulation of the initial guess.

My BC for sdot[0] and sdot[8]: reactor[0].v & reactor[1].v at the length 0 is reactor[0].v and reactor.v at the start &
reactor[0].p_ges & reactor[1].p_ges at the Point L (end) is eval to p_ges(at the end). But I dont know how I can programming it.

My BC for sdot[1:7] & sdot[9:15]: reactor[0].Fmol[i] at the start is the Inlet reactor[0].Fmol[i], same for reactor[1].Fmol[i ] & at the and is the change after the length is zero.

My BC for sdot[7] & sdot[15]: reactor[0].T and reactor[1].T at the Inlet is eval to the Inlet temperature & at the and is the change after the length is zero again like above.

Here is my Code:

from sub import substance, mixture
import numpy as np
from geo import reactor_volume, catalyst, membrane
from kinetics import kinetics
from math import pi
import matplotlib.pyplot as plt
import termplotlib as tpl
from scipy.integrate import solve_bvp as solve
from controller import controller

### Controller Settings
dpio = True
adiabatio = False
thermtubeio = True
thermsweepio = True
thermmemio = True
drainio = True
reacio = False

### Controller init
ctl = controller(dpio,adiabatio,thermtubeio,thermsweepio,thermmemio,drainio,reacio)
ctl.setcontroll()

## fixed variables
### Tube
Tref = 273 # reference temperature [K]
T0_t = 493 # temperature of tube inlet [K]
p0_t = 50*1e5 # pressure of tube inlet [Pa]
tube_diam = 0.016 # diameter of tube [m]
mdot_feed_t = 2.8*1e-5 # mass flow of tube inlet [kg/s]
reactor_length = 0.15 # length of reactor [m]

### Sweep
Rp = 0.1 # factor for controlling pressure in sweep
Rf = 0.1 # factor for controlling mass flow in sweep
T0_s = 400 # temperature of sweep inlet [K]
p0_s = Rp * p0_t # pressure of sweep [Pa]
outer_diam = 0.008 # outer diam only sweep
sweep_diam = tube_diam + outer_diam # diameter of sweep + tube [m]
mdot_feed_s = Rf * mdot_feed_t # mass flow of sweep inlet [kg/s]

###  Catalyst
material = "CuZnO"
kat_mass = 0.0348 # mass of catalyst [kg]
kat_diam = 0.0005 # diameter of catalyst particle [m]
kat_rho_part = 1775 # density of catalyst [kg/m3]

### Substances
Q = np.array([0, 0, 0, 0, 10**(-7), 0])
M = np.array([2e-3, 28e-3, 44e-3, 32.04e-3, 18e-3, 14e-3])
y0_t = np.array([0.82, 0.04, 0.03, 0, 0, 0.11])
y0_s = np.array([1, 0.0, 0.0, 0, 0, 0.0])
Pc = [20.5, 35*1e5, 73.83*1e5, 80.9464*1e5, 220.6*1e5, 34*1e5]
Tc = [43.6, 132.92, 304.21, 512.6, 647.096, 126.2]
Tb = [79, 81.7, 194.7, 337.668, 373.2, 77.355]
acent=[-0.216, 0.048, 0.2236, 0.556, 0.038, 0.038]
name = ["H2", "CO", "CO2", "MeOH", "H2O", "N2"]
pol = [False,True,False,True,True,False]
lambda_wall = 30

### object init 
sub_t, sub_s = [], []
rs = ["t","s"]

for j in range(len(rs)):
    for i in range(len(M)):
        exec(f"substances{j}{i}=substance(name[{i}],{i},pol[{i}],M[{i}],y0_{rs[j]}[{i}],Q[{i}],Pc=Pc[{i}],Tc=Tc[{i}],acent=acent[{i}],Tb=Tb[{i}])")
        exec(f"sub_{rs[j]}.append(substances{j}{i})") 

mix1,mix2 = mixture("tube", sub_t), mixture("sweep", sub_s)
gases = [mix1, mix2]
catalyst = catalyst(material,kat_diam,kat_mass,kat_rho_part)
kin = kinetics("Methanol", "Fromment")
tube, sweep = reactor_volume("tube",tube_diam,T0_t,p0_t,mdot_feed_t, mix1,catalyst,kin), reactor_volume("sweep",sweep_diam,T0_s,p0_s, mdot_feed_s,mix2)
reactor = [tube,sweep]
membrane = membrane(0,Q)

### Starting conditions
[exec(f"gases[{j}].setMolFrac(y0_{rs[j]})") for j in range(len(rs))]
[exec(f"gases[{j}].calcRho(T0_{rs[j]},p0_{rs[j]})") for j in range(len(rs))]
[exec(f"reactor[{j}].setMdot(mdot_feed_{rs[j]})") for j in range(len(rs))]
[exec(f"reactor[{j}].calcFmol(p0_{rs[j]},T0_{rs[j]})") for j in range(len(rs))]

for j in range(len(rs)):
    reactor[j].setFmol(reactor[j].Fmol_ges * reactor[j].gas.y)

### Solve
def xsolve(t,s):
    
    reactor[0].setPressure(s[0])
    reactor[0].setFmol(s[1:7])
    reactor[0].setTemp(s[7])
    
    reactor[1].setPressure(s[8])
    reactor[1].setFmol(s[9:15])
    reactor[1].setTemp(s[15])    
    
    sdot = [0 for _ in range(len(s))]
    
    for j in range(len(reactor)):
        for i in range(len(M)):
            reactor[j].gas.y[i] = reactor[j].Fmol[i]/reactor[j].Fmol_ges
    
    for i in range(len(M)):
        reactor[0].gas.substances[i].calcCpx(s[7])
        reactor[1].gas.substances[i].calcCpx(s[15])
        reactor[0].gas.substances[i].calcViscox(s[7],s[0])
        reactor[1].gas.substances[i].calcViscox(s[15],s[8])
        reactor[0].gas.substances[i].calcLambda(s[7])
        reactor[1].gas.substances[i].calcLambda(s[15])
    reactor[0].gas.calcRho(s[7],s[0])
    reactor[1].gas.calcRho(s[15],s[8])
    
    for i in range(len(reactor)):
        reactor[i].gas.calcCpm()
        reactor[i].gas.calcVisc()
    reactor[0].gas.calcLambdam(s[7])
    reactor[1].gas.calcLambdam(s[15])
    reactor[0].calcVelo(s[7],s[0],inner_diam=0)
    reactor[1].calcVelo(s[15],s[8],inner_diam=tube_diam)
    reactor[0].kat.calcPorosity(tube_diam)
    reactor[0].kat.calcRhobed()
    reactor[0].calcLength()
    reactor[1].length = reactor[0].length
    
    for i in range(len(reactor)):
        reactor[i].calcRe(kat_diam=kat_diam,eps=reactor[0].kat.porosity,inner_diam=tube_diam)
        reactor[i].calcPr()
        reactor[i].calcNu(inner_diam=tube_diam)
        reactor[i].calcAlpha_wall()
        
    reactor[0].kin.calcKeq(s[7])
    reactor[0].kin.calcHr(s[7],[reactor[0].gas.substances[i].cpx for i in range(len(M))])
    reactor[0].kin.calcReacrate(s[0],reactor[0].gas.y)
    
    membrane.setInterface(reactor[1].gas.y,reactor[0].gas.y,s[8],s[0],reactor[1].alpha_wall,reactor[0].alpha_wall)
    membrane.calcJ()
    membrane.calcU(sweep_diam,tube_diam,lambda_wall)
    
    sdot[0] = ctl.dpio * (-150 * reactor[0].gas.visc * (1 - reactor[0].kat.porosity)**2 * reactor[0].v / reactor[0].kat.porosity**3 / kat_diam**2 - 1.75 * (1 - reactor[0].kat.porosity) * reactor[0].gas.rho_T * reactor[0].v**2 / reactor[0].kat.porosity**3 / kat_diam)

    sdot[8] = 0
    
    for i in range(len(M)):
        sdot[i+1] = ctl.reacio * reactor[0].kat.rho_bed * np.pi * tube_diam**2 * reactor[0].kin.ri[i] - membrane.J[i] * 2 * np.pi * tube_diam
        sdot[i+9] = membrane.J[i] * 2 * np.pi * tube_diam
    
    sumJicpi = ctl.drainio * sum(membrane.J[i] * reactor[0].gas.substances[i].cpx for i in range(len(membrane.J)))
    sumFmcpitube = sum(reactor[0].Fmol[i] * reactor[0].gas.substances[i].cpx for i in range(len(reactor[0].Fmol)))
    sumFmcpisweep = sum(reactor[1].Fmol[i] * reactor[1].gas.substances[i].cpx for i in range(len(reactor[1].Fmol)))
           
    sdot[7] = ctl.thermtubeio * 2 * np.pi * tube_diam * (ctl.thermmemio * -membrane.memU * (reactor[0].T - reactor[1].T) - ctl.reacio* 0.5 * tube_diam* reactor[0].kat.rho_bed * (reactor[0].kin.dHr[0] * reactor[0].kin.reac_rate[0] + reactor[0].kin.dHr[1] * reactor[0].kin.reac_rate[1] + reactor[0].kin.dHr[2] * reactor[0].kin.reac_rate[2]) - sumJicpi * (reactor[0].T - Tref) ) / sumFmcpitube
        
    sdot[15] = ctl.thermsweepio * 2 * np.pi * tube_diam * (ctl.thermmemio * membrane.memU * (reactor[0].T - reactor[1].T) + sumJicpi * (reactor[0].T - Tref) ) / sumFmcpisweep
        
    return np.vstack(sdot)

def bc(pt,Fmolt[0],Fmolt[1],Fmolt[2],Fmolt[3],Fmolt[4],Fmolt[5],Tt,ps,Fmols[0],Fmols[1],Fmols[2],Fmols[3],Fmols[4],Fmols[5],Ts):

    return (np.array(), 
            np.array())

t_span = np.array([0, reactor_length])
t = np.linspace(t_span[0],t_span[1],50)
s0 = np.zeros((2, t.size))
s0 = np.array([])

sol = solve(xsolve,bc,t,s0,tol=1e-9,verbose=2)

If anyone can help me, I would be very grateful. I could also make the imported models available if necessary. Thank you very much for your attention.

Best regards

Understand and use the solve_bvp solver correctly