Hello to all,

I have been trying to model a simple diffusion problem consisting from 3 different domains that are connected through interfaces. Problem is just in 1D and and the domains are connected to each other in following manner.

```
domains = [x1 ∈ Interval(0.0, b)
x2 ∈ Interval(b,0.075)
x3 ∈ Interval(0.075,x_max)]
```

, where b=0.025 and x_max=0.275.

I have 3 equation system, each for one domain.

```
eq = [-D_1*Dxx1(u(x1)) + D_1*Vh/(8.314e3*T)*sigmaf_d_1(x1)*Dx1(u(x1)) + D_1*Vh/(8.314e3*T)*sigmaf_dd_1(x1)*u(x1) ~ 0
-D_2*Dxx2(v(x2)) + D_2*Vh/(8.314e3*T)*sigmaf_d_2(x2)*Dx2(v(x2)) + D_2*Vh/(8.314e3*T)*sigmaf_dd_2(x2)*v(x2) ~ 0
-D_1*Dxx3(h(x3)) + D_1*Vh/(8.314e3*T)*sigmaf_d_3(x3)*Dx3(h(x3)) + D_1*Vh/(8.314e3*T)*sigmaf_dd_3(x3)*h(x3) ~ 0]
```

Below the definition of the pdesystem and MOLFiniteDifference

```
@named pdesys = PDESystem(eq, bcs, domains, [x1, x2, x3], [u(x1), v(x2), h(x3)])
order = 2
discretization = MOLFiniteDifference([x1 => dx, x2 => dx, x3 => dx], nothing, approx_order=order)
prob = discretize(pdesys, discretization)
```

At this point the boundary conditions between interfaces are just limited to same concentration values. I’ll be implementing boundary energy related segregation condition later once this is solved.

Already after the first few runs the third domain work in a weird manner. The solution was completely linear although there exists relatively high pressure gradient and its second derivative which in this case are brought from readily interpolated values.

Why I started to think that not all the domains are considered is that once the discretization of the pdesystem is done, it returns only 2 element u0 vector, not 3 as I should expect in this case. Shown below

```
NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 2-element Vector{Float64}:
1.0
1.0
```

So my question is that, is the method of lines discretization limited to 2 domains only? I have to go for a different approach if that’s the case.

Thanks for the answers already in advance!