Finite difference Neumann boundary conditions: uneven weighting of edge nodes?

Originally asked this on math.stackexchange, but I figure it's also appropriate here. I'm reading through some finite difference code for a diffusion equation and came across something odd for the boundary conditions, and I was wondering about the validity of the method. I'm familiar with the finite difference method, but I've never seen this certain method before.

The equation I am solving is

$$\frac{\partial f}{\partial t}=C(x)\frac{\partial^2 f}{\partial x^2}$$

with boundary condition

$$\frac{\partial f}{\partial x}=0\ \text{on}\ x=0,1$$

across time, $t$, and position, $x$, with position-dependent diffusion constant, $C(x)$. It's being solved in an explicit way using the forward difference for time and the second-order central difference for space

$$\frac{f^\text{new}_i-f^\text{old}_i}{\Delta t}=C_i\frac{f^\text{old}_{i-1}-2f^\text{old}_{i}+f^\text{old}_{i+1}}{(\Delta x)^2}$$

along discretised space $\{x_1,\dots,x_n\}$. Rearranging the above equation gives

\begin{align} f^\text{new}_i=&f^\text{old}_i + \frac{C_i\Delta t}{(\Delta x)^2}\Big(f^\text{old}_{i-1}-2f^\text{old}_{i}+f^\text{old}_{i+1}\Big)\\ =&f^\text{old}_{i-1}(S_i)+f^\text{old}_{i}(1-2S_i)+f^\text{old}_{i+1}(S_i) \end{align}

where $S_i=\frac{C_i\Delta t}{(\Delta x)^2}$.

This is simply represented in the code. For the boundary at $x_1$ (and at $x_n$, but I'll only talk about the first as they are similar), there are some strange things going on. Usually I would do a Taylor series expansion of function $f$ at points $x_2$ and $x_3$

$$f_2=f_1+\frac{\partial f}{\partial x}\bigg\rvert_i \Delta x + \frac{1}{2} \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i (\Delta x)^2 + \mathcal{O}((\Delta x)^3)$$ $$f_3=f_1+\frac{\partial f}{\partial x}\bigg\rvert_i (2\Delta x) + \frac{1}{2} \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i (2\Delta x)^2 + \mathcal{O}((\Delta x)^3)$$

which, due to boundary condition $\frac{\partial f}{\partial x}=0$ can be rearranged to

$$\frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{2}{(\Delta x)^2}\Big(f_2 - f_1\Big)\tag{A}$$ $$\frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{1}{2(\Delta x)^2}\Big(f_3 - f_1\Big)\tag{B}.$$

Usually, I'd average $A$ and $B$ to determine how the boundary conditions are applied in the eventual matrix solution

$$\frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{1}{(\Delta x)^2}\Big(-\frac{5}{4}f_1 + f_2 + \frac{1}{4}f_3 \Big).$$

However, the code that I'm reading does something different and takes the boundary as a weighted sum of $A/5+4B/5$, as shown below

$$\frac{A+4B}{5} \Rightarrow \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{1}{(\Delta x)^2}\Big(-\frac{4}{5}f_1 + \frac{2}{5}f_2 + \frac{2}{5}f_3 \Big).$$

I don't think this as being inherently wrong, as the LHS is still $\frac{\partial^2 f}{\partial x^2}\big\rvert_i$, but I'm not sure. Is this a valid step to make? Is it ok to "weight" results from different nodes, and if so, in what situations would weighting or not weighting be useful?