Time-accurate simulations of the thermal interaction of flows and structures, also referred to as conjugate heat transfer (CHT), can be computationally expensive. Furthermore, given the multi-physics nature of many engineering problems, resolution of other coupled phenomena, in addition to CHT, may also be of interest. This thesis aimed at developing a flexible and efficient numerical procedure for solving unsteady (transient) conjugate heat transfer. High order time integration schemes are considered, in place of commonly used second order implicit schemes, to reduce the computational work of advancing the coupled problem in time. For flexibility, the partitioned method is adopted for solving the coupled problem. For strongly coupled problems, a strongly-coupled solution algorithm is presented where high order explicit first stage singly diagonally implicit Runge-Kutta (ESDIRK) schemes are used for time integration. For Dirichlet-Neumann conditions at the interface, stability and rate of convergence of subiterations at each stage are analyzed analytically. Based on the analysis, the domain with the higher effusivity is assigned the Neumann condition and the one with the lower effusivity the Dirichlet condition. Furthermore, the interface iterations converge with a rate approximately given by the ratio of thermal effusivities of the subdomains. For weakly coupled problems, an order preserving loosely coupled solution algorithm is presented in which a family of high order implicit-explicit (IMEX) Runge-Kutta schemes are used for time integration. The IMEX schemes consist of the ESDIRK schemes for advancing the solution in time within each subdomain, and equal order and number of stages explicit Runge-Kutta (ERK) schemes for explicit integration of part of the coupling terms. Based on a stability investigation, when the ratio of thermal effusivities of the subdomains is much smaller than unity, it is possible to take large Fourier numbers using the loosely coupled algorithm. Another topic studied in this thesis is the application of high order ESDIRK schemes to cell-centered collocated finite volume discretization of unsteady incompressible Navier-Stokes equations. In particular, a face-velocity interpolation procedure (Rhie-Chow) which preserves the temporal design order of the multi-stage ESDIRK schemes is introduced. The influence of iterative errors on temporal order is minimized by using an iterative time advancing algorithm.