Contact resistances play an important role in computational modelling of a wide range of physical phenomena, including heat transfer, electrical conductance, and electrophysiology. In all these cases, the systems considered consist of two or more sub-domains separated by interfaces, e.g. insulation layers, imperfect connections between adjacent surfaces, or cell membranes. Here we consider heat diffusion through a domain Ω = Ω_1 U Ω_2, where the two sub-domains are separated by an interface γ. The strong form of this problem is (1a) -div(k_i grad(u)) = f_i in Ω_i i=1,2 (1b) [k grad(u)] = 0, k grad(u) = R(u_1|_γ ,u_2|_γ)^(-1)[u] on γ. where the known, possibly non-linear function R describes the thermal resistivity of γ, and [.] is the jump operator. Ben Belgacem et al. [1] have previously considered a similar problem formulation, but restricted their analysis in two important ways: 1) The meshes of Ω_1 and Ω_2 must match on γ, 2) For 3D problems, R must be such that R^(-1) = a+b[u]^r with r=2 and a,b constants. For our application -- modelling multilayer insulation for use on large-scale liquid hydrogen carrier ships -- it is necessary to relax these assumptions. In a preliminary investigation, we have taken inspiration from Trujillo et al. [2] and Zheng et al. [3] to develop an iterative solution method for (1) that allows independent meshing of Ω_1 and Ω_2. Using the method of manufactured solutions and the radiation-like condition k grad(u) = u_1^4|_γ - u_2^4|_γ on γ, we have observed that our method converges in a series of numerical methods on both matching and non-matching meshes. Moreover, in a grid refinement study for matching meshes, we observe optimal convergence in the L2-norm for first-order local basis functions.