I have a partial differential equation of the form $$ \frac{d}{dt}f(x,t) = \Theta(x) f(x,t) \qquad \Theta(x) \sim \left[\frac{d^2}{dx^2} + k^2(x)\right] $$ subject to $f(x,t=0) = f_0(x)$, and $f(x=0,t) = f(x=1,t) = 0$. The function $k(x)$ is piecewise constant over $x\in[0,1]$. Call the $x$ intervals where $k$ is constant $I_j$. I would like to locally represent the solution $f$ as a sum of (say) Chebyshev polynomials within each interval $I_j$ and enforce continuity conditions at the boundaries of the intervals (continuity of $f$ and $\frac{\partial f}{\partial x}$).
Now, since the RHS operator is time independent, I could just solve this using an exponential: $f=e^{\Theta t}f_0$. My question is how I would compute the matrix elements of $\Theta$ to take into account the subdomain continuity conditions.
Alternatively, I could use a time marching method and make this a true domain decomposition method, but I still don't know how I would enforce the continuity conditions. I don't want to use a penalty method to enforce the constraints only approximately; I want them enforced as exactly as possible, and I really don't care about parallelism (so perhaps calling this DD is misleading). The literature on this probably exists, but I simply have no idea what I'm looking for.
I only give a rough description of $\Theta$, but it comes from the Helmholtz equation so it's hyperbolic, and a $2\times2$ block operator, with $f$ actually being a 2-vector of a function with its dual. This isn't terribly relevant to the question, since I just want to know how to enforce internal consistency.
