SummaryBy exploiting the meshless property of kernel‐based collocation methods, we propose a fully automatic numerical recipe for solving interpolation/regression and boundary value problems adaptively. The proposed algorithm is built upon a least squares collocation formulation on some quasi‐random point sets with low discrepancy. A novel strategy is proposed to ensure that the fill distances of data points in the domain and on the boundary are in the same order of magnitude. To circumvent the potential problem of ill‐conditioning due to extremely small separation distance in the point sets, we add an extra dimension to the data points for generating shape parameters such that nearby kernels are of distinctive shape. This effectively eliminates the needs of shape parameter identification. Resulting linear systems were then solved by a greedy trial space algorithm to improve the robustness of the algorithm. Numerical examples are provided to demonstrate the efficiency and accuracy of the proposed methods.