Why are models based on coded (-1,1 etc) and non-coded variables (as they are) very different? What should I use for publications? I am doing factorial experiments in R. I noticed when I use my variables as they are vs. coding them into -1,1, they are all very different. Here is my sample code. This dataset is adapted from the Experimentation for Improvement online course in Coursera.
#Example 1: Non-coded, as they are
runner1 <- data.frame(
  A = c(5, 10, 5, 10, 5, 10, 5, 10), 
  B = c(4, 4, 9, 9, 4, 4, 9, 9), 
  C = c(400, 400, 400, 400, 500, 500, 500, 500), 
  Y= c(39, 35, 40, 41, 40, 38, 41, 42))
runner1
runner1lm <- lm(runner1$Y ~ runner1$A*runner1$B*runner1$C)
runner1lm

#Example 2:Coded (-1, 1)

runner2 <- data.frame(
  A = c(-1, 1, -1, 1, -1, 1, -1, 1), 
  B = c(-1, -1, 1, 1, -1, -1, 1, 1),
  C = c(-1, -1, -1, -1, 1, 1, 1, 1), 
  Y= c(39, 35, 40, 41, 40, 38, 41, 42))
runner2
runner2lm <- lm(runner2$Y ~ runner2$A*runner2$B*runner2$C)
runner2lm

Which should I use? The coded or non-coded one? Also, how do I transform my entire dataset from coded to non-coded, if ever? Do I need a package for that? I was thinking about inputting coding parameters (like all 5 in column A will be -1, all 10 will be 1, so on and so forth like that) then letting R do its job and, boom!, I have a transformed data frame full of coded variables (except the dependent variable Y ofc). Basically I am asking if there's any code or package that can magically turn runner1 to runner2?

The models are identical.
One way to see this (although it's not a proof) is to compare their predictions and note that they are equal.
rbind(`First model`=predict(runner1lm), `Second model`=predict(runner2lm))


              1  2  3  4  5  6  7  8
First model  39 35 40 41 40 38 41 42
Second model 39 35 40 41 40 38 41 42


Of course the coefficients differ because you have used different methods to express the data with numbers.
We would expect the outputs to be interconvertible, though.  You can figure out how by writing the model predictions as formulas.  Let $A,B,C$ be the $\pm 1$ codes used in the second model.  In this model, expanding the interactions gives
$$y = \beta_0 + \beta_A A + \beta_B B + \beta_C C + \beta_{AB} A\times B + \beta_{AC}A\times C + \beta_{BC} B\times C + \beta_{ABC}A\times B \times C.$$
There are eight coefficients to estimate.  Their estimates, in the order given, are
$$(\beta_0, \beta_A, \ldots, \beta_{ABC}) = \left(39\frac{1}{2}, -\frac{1}{2}, 1\frac{1}{2}, \frac{3}{4}, 1, \frac{1}{4}, -\frac{1}{4}, -\frac{1}{4}\right),$$
as you saw when you printed runner2lm at the end of your code.
For instance, in the first record where $(A,B,C)=(-1,-1,-1),$ the prediction is
$$y = 39\frac{1}{2} -\frac{1}{2}(-1) + \cdots + -\frac{1}{4}(-1)(-1)(-1) = 39.$$
The first model changes how $A,$ $B,$ and $C$ are expressed.  Let's call these variables $a,b,$ and $c.$  The relationships are
$$\left\{\begin{aligned}
a &= \frac{15}{2} + \frac{5}{2}A;&\quad A &= (2a - 15)/5\\
b &= \frac{13}{2} + \frac{5}{2}B;&\quad B &= (2b - 13)/5\\
c &= 450 + 50C &\quad C &= (c - 450)/50.
\end{aligned}\right.$$
Plugging these expressions for $A,B,C$ into the model gives
$$y = \beta_0 + \beta_A(2a-15)/5 + \cdots + \beta_{ABC}\left((2a-15)/5\right)\left((2b-13)/5\right)\left((c-450)/50\right).$$
After expanding all these products and doing quite a bit of algebra, you find the corresponding coefficients for the first model.  For instance, the coefficient of $abc$ is $$\beta_{ABC}(2/5)(2/5)(1/50) = -4/(4\times 5\times 5\times 50) = 1/1250 = 0.0008.$$
You can inspect this directly:
tail(coefficients(runner1lm), 1)


runner1$A:runner1$B:runner1$C 
                       -8e-04


-8e-04 is computerese for $8\times 10^{-4} = 0.0008,$ agreeing with the calculation.
It is a worthwhile exercise to check some of the other coefficients in the same way: inspect the coefficient in one model; do the algebra to work out what it must be in the other model; and then confirm that through inspection.  Repeat this exercise until you feel you really understand the relationships between these two forms of the same model.
If you get stuck, simplify your example.  Start with a model of the form $y = \beta_0 + \beta_A A,$ with one explanatory variable.  Then move to a model with two variables and an interaction term.  After that you should have the hang of it.
Moral

When you report model estimates, make sure to explain how you chose numbers to express all the variables involved.

Some methods of numerical expression ("coding") lend themselves better to interpretation than others, but any method that creates a mathematically predictable one-to-one correspondence, as examined above, will work.
Computing bonus
Because the relationships between these forms of numerical expression can be computed algebraically, they permit easy, efficient transformation between the data frames.  For instance, the expression
15/2 + (5/2)*runner2$A

converts the A variable in runner2 (which I named $A$ in the formulas above) to the A variable in runner1 (which I named $a$).  Thus, you don't need to create any translation tables or use any special packages: addition and multiplication take care of everything.
