Recursive Implementation of the Gaussian Filter
The answer was simple, the article uses the coefficients value on one hand where the MATLAB implementation on the other.
Namely, a minus sign should be added.

Here's the correct code:

qFactor = 5;

b0Coeff = 1.57825 + (2.44413 * qFactor) + (1.4281 * qFactor * qFactor) + (0.422205 * qFactor * qFactor * qFactor);
b1Coeff = (2.44413 * qFactor) + (2.85619 * qFactor * qFactor) + (1.26661 * qFactor * qFactor * qFactor);
b2Coeff = (-1.4281 * qFactor * qFactor) + (-1.26661 * qFactor * qFactor * qFactor);
b3Coeff = 0.422205 * qFactor * qFactor * qFactor;

normalizationCoeff = 1 - ((b1Coeff + b2Coeff + b3Coeff) / b0Coeff);

vDenCoeff = [b0Coeff, -b1Coeff, -b2Coeff, -b3Coeff] / b0Coeff;


vXSignal = zeros(61, 1);
vXSignal(31) = 10;

vYSignal = filter(normalizationCoeff, vDenCoeff, vXSignal);
vYSignal = filter(normalizationCoeff, vDenCoeff, vYSignal(end:-1:1));

figure();
plot(vYSignal);