It is easy to show that if $f\colon[0,1]\to\mathbb R$ and $|f|\leq A$ and $|f''|\leq B$ then~$|f'|\leq 4A+B$. Indeed, by Taylor formula with remainder $f(x)=f(c)+(x-c)f'(c)+\frac12(x-c)^2f''(d)$ where $d$ is between $x$ and $c$. Therefore $f'(c)=\frac{f(x)-f(c)-\frac12(x-c)^2f''(d)}{x-c}$. Now let $c$ be a maximum of $|f'|$. We can assume without loss of generality that $c\geq\frac12$ (otherwise apply the same argument to the function $f(\frac12-x)$ to get the same bound). Then \begin{aligned} |f'(c)|&\leq2\left|f(0)-f(c)-\frac12(0-c)^2f''(d)\right| \\&\leq2|f(0)|+2|f(c)|+c^2 |f''(d)| \\&\leq4A+B \end{aligned} as required. In particular, $f$ is $(4A+B)$-Lipschitz.
The question is, can one obtain similar bounds for the Lipschitz constant of $f$ if one does not require differentiability but only assumes a bound on second differences: $f(x)-2f(x+h)+f(x+2h)\leq Bh^2$ ?
Yes, the supremum $A = \sup_x |f(x)|$ and the second divided difference $$B = \sup_{x,h} \left|\frac{f(x-h)-2f(x)+f(x+h)}{h^2}\right|$$ control the Lipschitz constant of $f$. One way to see this is to use mollification: convolve $f$ with a bump function $\phi_\epsilon$ supported on $[-\epsilon,\epsilon]$ (so the convolution will be defined on $[\epsilon,1-\epsilon]$), and observe that both bounds are inherited by the mollified function $f_\epsilon = f*\phi_\epsilon$. Since $f_\epsilon$ is smooth, the proof you have works for it. As $\epsilon\to 0$, we have $f_\epsilon\to f$ in $L^1$, and the set of Lipschitz functions with uniformly bounded Lipschitz constant is closed under $L^1$ convergence.
Alternative approach, without convolution
Suppose that for some $a, h$ the ratio $|f(a+h)-f(a)|/h$ is large (say, $=L$). WLOG assume $(f(a+h)-f(a))/h>0$. The bound on second differences tells us that $$ \left| \frac{f(a+2h)-f(a+h)}{h} - \frac{f(a+h)-f(a)}{h} \right| \le Bh $$ So, the divided differences $(f(a+kh+h)-f(a+kh))/h$ will remain over $L/2$ while $k<L/(2Bh)$. This gives us $N = L/(2Bh)$ consecutive subintervals* of size $h$, on each of which $f$ grows by more than $Lh/2$. So the total growth of $f$ is over $NLh/2 = L^2/(4B)$, which implies $L^2/(4B)\le 2A$, hence $L\le \sqrt{8AB}$.
(*) Boundary effect: we may run out of the interval $[0,1]$ in the process. The geometric-mean bound $L\le \sqrt{8AB}$ is valid for functions on the entire line (or half-line), but not on a bounded interval. Within the interval $[0,1]$ we are limited to about $1/h$ steps. If $L/(2Bh) > 1/h$, we end up with $1/h$ subintervals, hence the growth of $f$ is $L/2$, and the result is $L/2 \le 2A$, hence $L\le 4A$. Final outcome: $$ L\le \max(4A, \sqrt{8AB}) \tag{1} $$
Note on rounding. One may worry about rounding the quantities like $1/h$ one way or another, and the effect it has on the estimate (1). But notice that $h$ can be taken arbitrarily small by subdividing the initial subinterval $[a,a+h]$, so rounding does not matter.