Radial basis functions (RBF) are widely used in many areas especially for interpolation and approximation of scattered data, solution of ordinary and partial differential equations, etc. The RBF methods belong to meshless methods, which do not require tessellation of the data domain, i.e. using Delaunay triangulation, in general. The RBF meshless methods are independent of a dimensionality of the problem solved and they mostly lead to a solution of a linear system of equations. Generally, the approximation is formed using the principle of unity as a sum of weighed RBFs. These two classes of RBFs: global and local, mostly having a shape parameter determining the RBF behavior. In this contribution, we present preliminary results of the estimation of a vector of �optimal� shape parameters, which are different for each RBF used in the final formula for RBF approximation. The preliminary experimental results proved, that there are many local optima and if an iteration process is to be used, no guaranteed global optima are obtained. Therefore, an iterative process, e.g. used in partial differential equation solutions, might find a local optimum, which can be far from the global optima. © Springer Nature Switzerland AG 2020.