This work focuses on interface problems that arise in many areas of science and engineering. In that problems some parameters of the model exhibit jump in their values at the surface (3D), curve (2D) or point (1D) which is called interface. A mathematical model of such problems is in the form of differential equation completed by the boundary conditions and the conditions specified on that interface. The discontinuities of the model parameters may lead to strong or weak discontinuities of the solution across the interface. In this work the domain of the problem is divided by the interface or interfaces into subdomains, on which the mathematical model of the interface problem is defined. The model is originally defined in the local form and then is transformed to the global multicriterial optimization problem for the whole domain in the form of the convex functional to be minimized. The functional is composed of sum of squares of the residues of the governing equation, boundary and interface conditions, where the squares are expressed by the inner products. The components of the functional have different physical dimensions, but the physical consistency of the functional is affirmed. The minimization process leads to a variational problem, which in the finite space has the form of system of algebraic equations with the degrees of freedom as unknowns. The discrete problem is symmetric and well-defined, even for high-order approximation. In the discrete form of the interface problem the inner products are approximated by their discrete counterparts obtained by the collocation procedure. The usage of the collocation points together with the meshless approximation places this method among the truly meshless techniques, since no mesh structure is used for the approximation nor the numerical integration. The basis functions, which are applied for the approximate solution, are defined separately in each of the subdomains, and so different types basis functions can be applied in those subdomains. Although, the discrete problem is well-defined, it is a well-known property, that high-resolution approximation can spoil the condition of the matrix of the problem. To avoid such a problem the special null-space solver has been employed. To examine the method, several examples of interface problems have been solved showing its high rate of convergence and stability.