As the use of computer simulation for scientific discovery increases there is a growing need for reliable multiphysics simulations. Although an exact definition for multiphysics problems is difficult to state, these problems include simulations where two or more component models are coupled to simulate events beyond either individual component. Aware of the growing prevalence of multiphysics simulations, this work identifies potential impediments to efficient and stable computation and proposes procedures to address the concerns. We first introduce the implicitly coupled multiphysics framework, through which we limit the domain of problems we consider, and Jacobian-free NewtonKrylov methods, which will be the primary technique for solving these nonlinear systems. Kernel-based approximation theory is introduced, providing a method for coupling component models with different discretizations. These topics have both theoretical and computational implications which are later applied in the context of multiphysics simulations. Our first contribution is the analysis of preconditioning nonlinear systems produced during multiphysics simulations. The motivating application is computational plasma physics, specifically, the edge region of a tokamak reactor performing magnetic confinement fusion. For sufficiently difficult parameter choices, existing preconditioners have proved ineffective or inefficient at producing useful Newton steps. We analyze the various components of this simulation and determine why preconditioners acting on the entire simulation do not perform well. Using these insights we develop an operator-specific (or physics based) preconditioner which allows for better performance while improving parallel scalability. Preconditioning is an important component of the Jacobian-free NewtonKrylov method because it allows a faster solution of the Newton search directions. Equally critical, although less susceptible to poor performance, are the Jacobianvector products within the linear solvers which are approximated using Taylor series. We study this process in the multiphysics setting and propose improvements that allow for greater accuracy when two components on different scales are coupled together. The final issue we consider in multiphysics systems also receives the most treatment. Simulating a multiphysics system requires coupling between the individual components, and in this thesis we discuss the use of kernel-based scattered data interpolation to perform the coupling. A new technique, based on an orthonormal expansion of the kernel, is developed which allows us to evaluate meshfree radial basis interpolants in arbitrary dimensions without the ill-conditioning often present for accurate kernel choices. These eigenfunctions are derived using Hilbert-Schmidt theory, and tested on interpolation problems in up to five dimensions. After the eigenfunction method is validated, we show how it can be used to approximate derivatives of a function given only scattered data. Once this is proven, results are given for a multiphysics simulation using meshfree interpolation. These results are compared to the standard discretization scheme for interpolating between models, with higher order results possible using meshfree interpolation. Additionally, because the meshfree approach works with scattered data, it provides a more general method for coupling two models with mismatched grids. In addition to the viability of meshfree coupling for multiphysics, we also consider the computational cost. The preconditioning discussion from earlier is applied here to choose a good preconditioner for the fully coupled system, which is important for reducing the cost of linear solves. The remaining content focuses on other aspects of the eigenfunction expansion which are relevant to the wider computational community. Using the derivative approximation method created for coupling multiphysics components, we show how boundary value problems can be solved by collocation using the eigenfunction basis. We also explore the use of the eigenfunction technique in the Method of Particular Solutions, which demonstrates the benefit of solving elliptic problems with a joint collocation/Method of Fundamental Solutions approach. The final results deal with a statistical framework for determining appropriate parameterizations of kernels, which is necessary to realize the optimal interpolation accuracy discussed earlier. These methods have previously suffered from the ill-conditioning addressed by the eigenfunction approach, and, using the new stable basis, we reconsider their viability as predictive tools.