I would like to suggest a somewhat different approach compared to the other answers, though @barron has indirectly discussed the same thing.
Instead of optimizing your function directly, i.e. by evaluating it at a series of points $\mathbf{x}_1,\mathbf{x}_2,\ldots,\mathbf{x}_k$ points that (hopefully) converge to a (local) optimum, you could use the concept of $\textit{surrogate modelling}$, which is very well suited for problems of the type you describe (high cost, smooth, bounded, low dimensional i.e. less than 20 unknowns).
Specifically, surrogate modelling works by setting up a model function $c \in \mathbb{R}^d \rightarrow \mathbb{R}$ of your true function $f \in \mathbb{R}^d \rightarrow \mathbb{R}$. The key is that while $c$ of course does not perfectly represent $f$, it is far cheaper to evaluate.
So, a typical optimization process would be as follows:
- Evaluate $f$ at a set of $j$ initial points $\mathbf{x}_1,\mathbf{x}_2,\ldots,\mathbf{x}_j$. Note that derivatives are not needed. Also note that these points should be distributed evenly throughout the search space, e.g. by Latin Hypercube Sampling or a similar space-filling design.
- Based on this original dataset, create a model function $c$. You could use cross validation to validate your model (i.e. use only a subset of the original $j$ points to create $c$, and then use the remainder of the dataset to check how well $c$ predicts those values)
- Use a criterion such as the Expected Improvement (EI) criterion to find out where to ''fill in'' more samples to make $c$ more accurate by sampling $f$. This is actually far better studied theoretically than it might seem, and the EI criterion is very well researched. The EI criterion is also not a greedy criterion, so you both get good overall improvement of the model accuracy, whilst prioritizing accuracy near potential optima.
- If your model is not accurate enough, repeat step 3, else use your favourite optimization routine to find the optimum of $c$, which will be very cheap to evaluate (so you could use any routine you want, even ones that requires derivatives, or simply evaluate the function in a fine mesh).
In general, this is what is meant by EGO, Efficient Global Optimization, as @barron suggested. I would like to stress that for your application, this seems perfectly suitable — you get a surprisingly accurate model based on relatively few evaluations of $f$, and can then use any optimization algorithm you want. What's often also interesting is that you can now evaluate $c$ on a mesh and plot it, thereby gaining insight into the general appearance of $f$. Another interesting point is that most surrogate modelling techniques also provide statistical error estimates, thereby allowing uncertainty estimation.
How to construct $c$ is of course an open question, but often Kriging or so-called space-mapping models are used.
Of course, this is all quite a bit of coding work, but a lot of other people have done very good implementations. In Matlab, I only know of the DACE software toolbox DACE is free. TOMLAB might also offer a Matlab package, but costs money — however, I believe it also works in C++ and has far more capabilities than DACE will ever have. (Note: I am one of the developers of the new version of DACE, soon to be released, which will offer added support for EGO.)
Hope that this rough overview has helped you, please ask questions if there are points that can be made more clear or stuff I've missed, or if you would like further material on the subject.
