I am trying to solve the Heat equation with Robin Boundary condition:
$$ u_t(x,t) = u_{xx}(x,t), \\ u(x,0) = g(x), \\ u(0,t) + u_x(0,t) = h_0(t), \\ u(1,t) + u_x(1,t) = h_1(t)$$
for $ 0\leq x\leq1$ and $t>0$, by using twice the composition of a fourth-order compact finite difference (CFD) approximation of the first derivative (taken from S. Lele (1991)) for the spatial discretization, i.e.
$$ M \; \partial_x \textbf{u} = A \;\textbf{u} \\ M \; \partial^2_x \textbf{u} = A \; \partial_x \textbf{u}, $$
where $M, A$ are two sparse matrices, corresponding to the coefficients obtained from the CFD scheme. This would simplify as follows
$$ \partial^2_x \textbf{u} = (M^{-1} A)^2 \;\textbf{u}. $$
As adviced in S. Lele (1991), I use an upwinding and downwinding scheme at the boundaries, so I avoid using ghost cells with the derivatives.
I am aware this could be done in an easier manner by using an approximation of the second derivative, as done in J. Malele et al. (2022), but my problem is in fact more complicated than the one shown here and requires me to use twice the composition of the first derivative.
And here, my question: How could I include the Robin boundary conditions?
I was advised to include the equations as additional lines on the last system, so that
$$ \hat{\textbf{u}} = \begin{bmatrix} u_x\Big|_1 \\ u_0 \\ u_1 \\ \vdots \\ u_N \\ u_x\Big|_{N-1} \end{bmatrix} $$
would be my new unknown, however, that makes a new matrix appear on the left hand side of the equation, i.e. left-multiplying $ \partial_x^2 \hat{\textbf{u}} $, which is a singular matrix. The fact that this is a singular matrix would then complicate the resolution of the ODE with another method time integrator.
Thank you for your time.
P.S.: The matrix $M$ is the matrix of coefficients of the derivatives. That means that if the schemes are of this form (as per S. Lele (1991)):
$$ \alpha f_{i-1}' + f_i + \alpha f_{i+1}' = \frac{a}{2h} \left( f_{i+1} - f_{i+1} \right) \quad \forall i = 1,2,\dots, N-1$$
at the inner nodes and
$$ f_0' + \beta f_1' = \frac{1}{h}\left( bf_0 + cf_1 +df_2 + ef_3 \right) $$
at the left boundary, i.e. $i = 0$, and something similar for the right boundary, then the Matrix of coeffients will be:
$$ M := \begin{bmatrix} 1 & \beta & 0 & 0 & \dots \\ \alpha & 1 & \alpha & 0 & \ddots \\ 0 & \alpha & 1 & \alpha & \ddots \\ \vdots & \ddots & \ddots & \ddots & \ddots \end{bmatrix} $$