So I have managed to find two approaches for solving the coupled system:

1. Solve equations (2) - (4) using a standard ODE solver with an initial condition for $f_0=f(0,T)$ and $g_0=g(0)$. The newly calculated value of $q_0$ can then be substituted into equation (1) to find $f$. Use fixed point iteration with under-relaxation until the residuals for $f$ and $g$ converge and then progress to the next space step

2. Convert the system to a DAE system by moving everything to the RHS in (2-5) and discretising all terms in time using finite differences. This can be solved using a DAE solver, and specifying a mass matrix with zeros along the diagonal for (2-5).

Both methods work, however method (1) takes many iteration to converge, even with under-relaxation and method (2) becomes extremely demanding on memory when there are many time points

So I have managed to find two approaches for solving the coupled system:

1. Solve equations (2) - (4) using a standard ODE solver with an initial condition for $f_0=f(0,T)$ and $g_0=g(0)$. The newly calculated value of $q_0$ can then be substituted into equation (1) to find $f$. Use fixed point iteration with under-relaxation until the residuals for $f$ and $g$ converge and then progress to the next space step

2. Convert the system to a DAE system by moving everything to the RHS in (2-5) and discretising all terms in time using finite differences. This can be solved using a DAE solver, and specifying a mass matrix with zeros along the diagonal for (2-5).

Both methods work, however method (1) takes many iteration to converge, even with under-relaxation and method (2) becomes extremely demanding on memory when there are many time points

Following suggestions in the comments, another method is to discretise in $X$, rather than in $T$ and time-march, leading to the following ODE system (after some re-arrangement):

$$\frac{\partial f}{\partial T} = \gamma$$ $$\frac{\partial}{\partial T}(\frac{f^2}{2} + B \gamma) = \frac{\partial f}{\partial X} + ...$$ $$ ... $$

where $\gamma$ is a new variable to remove the 2nd order derivative in (1). However, I am unsure how to treat the term: $\frac{f^2}{2} + B \gamma$ on the LHS. Does anyone have any suggestions?