So Marshak boundary conditions are mixed

\left . \left ( E-\frac{2}{3\kappa}\nabla{E} \right ) \right |_0 = \frac{4}{c} F

One could instead do

E=\frac{4}{c}F

or

\frac{c}{3\kappa} \nabla E = F

If we discretize this equation we get (E_g is ghost zone, E_i is internal zone at left edge)

E_g + E_i - \frac{4}{3\kappa \Delta x} \left ( E_i - E_g \right ) = \frac{8}{c} F

which we can solve for E_g as

E_g = \frac{\left ( 1 - \frac{3\kappa \Delta x}{4} \right ) E_i + \frac{3\kappa \Delta x}{4} \left ( \frac{8}{c} F \right )}{1 + \frac{3 \kappa \Delta x}{4}}

and we have E_g - E_i = \frac{\left ( -\frac{6\kappa \Delta x}{4} \right ) E_i + \frac{3\kappa \Delta x}{4} \left ( \frac{8}{c} F \right )}{1 + \frac{3 \kappa \Delta x}{4}}

and

E_g + E_i = \frac{ 2 E_i + \frac{3\kappa \Delta x}{4} \left ( \frac{8}{c} F \right )}{1 + \frac{3 \kappa \Delta x}{4}}

If we then plug this into the time evolution equation:

E^{n+1}_i-E^{n}_i = \alpha \left ( E_g^*-E_i^* \right )

we get

E^{n+1}_i-E^{n}_i = \alpha \left ( \frac{-\frac{6\kappa \Delta x}{4} E_i^* + \frac{3\kappa \Delta x}{4} \left ( \frac{8}{c} F \right )}{1 + \frac{3 \kappa \Delta x}{4}} \right )

where E^* = \psi E^{n+1} + \bar{\psi} E^n is the time averaged energy and depends on the time stepping (crank-nicholson or forward euler)

and finally we arrive at

\left ( 1+\frac{2 \psi \alpha}{1+\frac{4}{3\kappa\Delta x}} \right ) E^{n+1}_i = \left ( 1-\frac{2 \bar{\psi} \alpha}{1+\frac{4}{3\kappa\Delta x}} \right ) E^{n}_i+ \frac{\frac{8\alpha}{c} F}{1+\frac{4}{3\kappa \Delta x}} This is very similar to what we normally get ie

\left ( 1+\psi \alpha \right ) E^{n+1}_i - \psi \alpha E^{n+1}_g = \left ( 1-\bar{\psi} \alpha \right ) E^{n}_i + \bar{\psi} \alpha E^{n}_g

and noting that \alpha = \frac{c \Delta t}{3 \kappa \Delta x^2} for the limiter \lambda=1/3

or flux is just

F \frac{\Delta t}{\Delta x} \rightarrow F \frac{\Delta t}{\Delta x} \frac{2}{1+\frac{3\kappa\Delta x}{4}}

Note that E^{n+1}_i = E^{n}_i for

E = \frac{4F}{c}