How to add Gauge/Tau Method

Im trying to use the interlace operator to get the US matrix representation of the system for functions $u,v, q$ and real numbers $\tau_1, \tau_2$
kl

$$
\begin{matrix}
u^{(4)} + q = f_1 \
v^{(4)} + q^{(1)} = f_2 \
-u^{(1)} + v + \tau_1 \phi_1 + \tau_2 \phi_2 = f_3 \
u(\pm 1) = u_{\pm} \
v(\pm 1) = v_{\pm}
\end{matrix}
$$

The thing is that I'm not sure how to add this variables $\tau_1, \tau_2$ through the block structure that comes in approx fun. If there was not any dummy variables, I was writing the system as

a,b = -1,1;
d = a..b;
D⁴ = Derivative(d, 4)
D¹ = Derivative(d, 1)

γ₀ = Dirichlet(d)
γ₁ = Neumann(d)
operator = [ γ₀  0    0 ;
             0   γ₀   0;
            γ₁  0    0 ;
             0   γ₁   0 ;
             I  -D¹  0;
            D⁴  0  I;
             0   D⁴  D¹
            ]

Another issue, is that there is a missmatch on how many coefficients I can enforce each of the equations, and not sure how to resolve that either using the notation of the library

Start by tracing the approx fun block structure and interlace operator APIs referenced in the issue. Work through how the system represents the dummy variables τ₁ and τ₂ and how many coefficients each equation can enforce. Done means a reproducible representation of the displayed system with the coefficient-count mismatch resolved.