In this paper we deal with a numerical method for a complex free boundary value problem (FBP) in 1D, arising from chemical kinetics and involving two moving internal boundary points and six unknown concentration profiles. Basically, the method consists of three steps: (1) a suitable fixed domain transformation for each of the three time-varying intervals, which results in a strongly nonlinear, nonlocal boundary value problem (BVP); (2) a nonstandard central difference method with respect to the space variable, that takes properly into account the various transition conditions, in particular by using quadratic Lagrange interpolation to approximate the involved concentration profiles near the relevant grid points; (3) a time discretisation method for the resulting initial value problem (IVP) for the nonlinear system of 1st order ODEs, that takes fully profit of the special structure of the mass matrix, implying that only a full \(8\times 8\) submatrix must be inverted at each discrete time point. A numerical example is presented.