So I have a system of coupled PDEs that represent the movement of an interface of fluid. The interfaces are functions of time and space, ie. $h_1(x,t)$, $h_2(x,t)$, $h_3(x,t)$. The PDEs are horrendous (so I won't write in full here) but they are of the form:

$h_{1_t} = f(h_1,h_2,h_3,h_{1_x},h_{2_x},h_{3_x},h_{1_{xx}},h_{2_{xx}},h_{3_{xx}}) \\ h_{2_t} = g(h_1,h_2,h_3,h_{1_x},h_{2_x},h_{3_x},h_{1_{xx}},h_{2_{xx}},h_{3_{xx}}) \\ h_{3_t} = p(h_1,h_2,h_3,h_{1_x},h_{2_x},h_{3_x},h_{1_{xx}},h_{2_{xx}},h_{3_{xx}})$

I wish to solve in a domain $|x| \leq 0.5$ and for time say $t \in (0,1)$. The initial conditions are as follows:

$h_1(x,0) = a(x) \quad \text{for} \, -0.5 < x \leq 0.1 \\ h_2(x,0) = b(x) \quad \text{for} \, 0 \leq x \leq 0.1 \\ h_3(x,0) = c(x) \quad \text{for} \, 0 \leq x < 0.5$

$a(x), b(x), c(x)$ are known functions and join up continuously, such that $a(0.1) = b(0.1)$ and $b(0) = c(0)$. Also $a(-1/2) = c(1/2) = 0$ to conform with the following boundary conditions:

$h_1(-0.5,t) = 0 \\ h_3(0.5,t) = 0 \\ h_2(\max(x),t) = h_3(\max(x),t) \\ h_2(\min(x),t) = h_1(\min(x),t)$

such that the ends of $h_1$ and $h_3$ are pinned at the end of my domain, and $h_2$ is always smoothly connecting $h_1$ and $h_3$.

I wish to solve this system numerically with NDSolve or NDSolveValue. However, I am having some issues with implementation. I am unsure as to how to input the discontinuous boundary conditions?

My current go is something like this:

NDSolve[{
   D[h1[x, t], t] == stuff1,
   D[h2[x, t], t] == stuff,
   D[h3[x, t], t] == stuff,
   h1[x,0] == a[x], h2[x,0] == b[x], h3[x,0] == c[x], h1[-0.5,t] == 0, h3[0.5,t] == 0
  },
 {h1,h2,h3}, {x,-0.5,0.5}, {t,0,1}
]

Clearly this isn't working, or handing the $h_2$ boundary condition. The functions $a(x), c(x), b(x)$ are all InterpolatingFunction running for their respective ranges (ie. $a(x)$ is defined from $-0.5 < x < 0.1$ etc).

Any ideas on how to deal with this type of boundary condition?