approximating a function without derivatives

I have a function of 3 variables $$F(a, b, c)$$ , and want to decompose the change in F into changes due to individual changes in $$a,b,c$$. Need to do it WITHOUT calculating actual derivatives - e.g. if the function is unknown but you have a blackbox that outputs $$F(a,b,c)$$ for inputs $$a,b,c$$.

At first, I thought of this: $$F(a,b,c) - F( a1, b1, c1) = [F(a,b,c) - F( a1, b, c)] + [F(a1,b,c) - F( a1, b1, c)] + [F(a1,b1,c) - F( a1, b1, c1)]$$

However, this doesn't really give a reasonable expectation of changes in F due to changes in a variable and it also depends on the order of change in variables.

So I am doing this:

$$F(a,b,c) - F( a1, b1, c1) \approx [F(a,b,c) - F( a1, b, c)] + [F(a,b,c) - F( a, b1, c)] + [F(a,b,c) - F( a, b, c1)]$$

I think this is equivalent to a taylor expansion where you ignore all the cross terms, but DON'T ignore the higher order terms of single variables.

Without any way of estimating the cross terms, is there a better way to do this? Or is there any way to estimate the approximation error term with simplest possible approximations but still without actual derivatives?

The reason to do this is to separate the changes in F due to variables that are more forecastable vs those that are completely random.

For illustration, let us start with two variables $$(a,b)$$ and you want to approximate the function by a plane. This means that, locally, the model is linear with respect to $$a$$ like $$F(a,b)= \alpha(b) + \beta(b) \,a$$ But, to be linear in the other dircetion means that $$\alpha(b)=\alpha_0+\alpha_1\,b \qquad \text{and}\qquad \beta(b)=\beta_0+\beta_1\,b$$

If you recombine everything, then the approximation is $$F(a,b)=A + B \,a+C\,b+D\,a\,b$$

If we extend to $$(a,b,c)$$ it will be $$F(a,b,c)=A + B \,a+C\,b+D\,c+E \,a\,b+F \,a\,c+G\,b\,c+H\,a\,b\,c$$ which takes into account the cross-interactions of the variables.