A Multi-Material Pin-Fin Heat Dissipation Modeling Method Based on an Efficient Shooting–Secant Hybrid Strategy and High-Order Numerical Integration

Wance Chen^*, Xinyue Zhang

Manchester Metropolitan Joint Institute, Hubei University, Hubei University, Wuhan, China, 430062

*Corresponding author: 202231123003011@stu.hubu.edu.cn

Keywords: Pin-Fin Heat Dissipation, High-Order Numerical Integration, Secant Method, Multi-Material Modeling

Abstract: With the continuous increase in power density of integrated circuits, microscale heat dissipation bottlenecks have become a critical factor limiting the performance of electronic systems. Pin-fin heat sinks, owing to their high surface-to-volume ratio and mature manufacturing processes, exhibit significant advantages in high heat flux density scenarios. However, the temperature distribution characteristics of pin fins under different material conditions, as well as the stability and accuracy of their numerical solutions, have not been systematically investigated. In this work, a unified numerical framework combining the shooting method with the secant method is proposed to address steady-state heat conduction boundary value problems (BVPs) by transforming them into initial value problems (IVPs). A derivative-free iterative scheme is employed to achieve convergent estimates of the initial derivative efficiently. For numerical integration, both the low-order Euler method and the fifth-order Runge-Kutta method based on Dormand-Prince coefficients are utilized to accurately reconstruct the temperature field of the IVP. Comparisons with analytical solutions demonstrate the high accuracy and robustness of the proposed method in multi-material systems. Simulation results reveal a strong positive correlation between the tip temperature and the material thermal conductivity. Further simulations on copper, aluminum, steel, and titanium validate the effectiveness of the method. This framework not only provides a quantitative computational basis for the thermal design of multi-material pin-fin heat sinks but also offers a generalizable algorithmic approach for the high-precision solution of complex boundary value problems.

1.Introduction

With the continuous rise in the power density of modern integrated circuits, thermal management has become a critical bottleneck limiting further performance enhancement of electronic systems. As transistor sizes continue to shrink and functional integration increases, local heat flux densities can reach hundreds of W/cm², resulting in significant hotspot effects and temperature non-uniformities. Traditional cooling methods mainly rely on increasing the convective heat transfer coefficient or enlarging the heat dissipation area. However, under microscale conditions, the limited heat transfer interface and thermal conductivity path make it difficult to meet the demands under high heat flux densities. This not only leads to reduced computational efficiency but also introduces reliability risks such as material fatigue and electromigration. Authoritative reviews have extensively summarized this issue. Wang et al. systematically reviewed the thermal bottleneck problem in three-dimensional stacked chips, pointing out that with the enhancement of hotspot heat flux and the limitation of heat dissipation paths, there is an urgent need for innovative and efficient heat dissipation strategies [1]. Garimella et al. highlighted that next-generation electronic systems face long-term challenges such as high heat flux densities, hotspot control, and multi-physics coupling [2]. These challenges significantly limit device performance and reliability. These findings collectively indicate that the development of compact and efficient heat dissipation structures holds significant academic and engineering value. Among various passive cooling technologies, pin-fin heat sinks have attracted widespread attention due to their high surface-to-volume ratio, mature manufacturing processes, and compatibility with existing packaging technologies. Their compact geometric structure effectively expands the heat transfer interface and enhances convective heat dissipation capabilities under limited space conditions. As a result, pin-fin structures are widely used in high-power-density scenarios such as microprocessors, power electronic devices, and optoelectronic systems. Existing research primarily focuses on optimizing the geometric parameters of pin fins, such as diameter, length, and arrangement, to reduce overall thermal resistance and minimize pressure drop. Numerical simulations and experimental studies have revealed the effects of these factors on heat transfer enhancement and flow characteristics. For example, computational fluid dynamics studies indicate significant differences in vortex formation and local heat transfer coefficients between staggered and inline arrangements, while experimental results validate the scalability of pin-fin modules in compact electronic devices. However, existing research still has limitations: firstly, most studies primarily focus on geometric structure or convective heat transfer, while the thermal conductivity characteristics of pin fins under different material conditions have not been systematically analyzed; secondly, in numerical modeling, low-order numerical integration methods often lack accuracy when solving the second-order differential equations of pin fin thermal conduction, especially when the boundary conditions are distributed at the ends of the domain. Therefore, there is an urgent need for a numerical framework that ensures stability, accuracy, and computational efficiency under multi-material and multi-operating conditions.

This study addresses the steady-state heat conduction boundary value problem for pin-fin structures and proposes a unified numerical framework for its solution. The method transforms the boundary value problem into an initial value problem, combining the secant method to achieve rapid convergence estimates of the unknown initial derivatives, ensuring iterative stability without the need for derivative information. For numerical integration of the temperature field, the framework simultaneously employs the low-order Euler method and the fifth-order Runge-Kutta method based on Dormand-Prince coefficients, balancing computational efficiency and accuracy. To further enhance the algorithm’s robustness, the step size and convergence threshold were optimized, and its high-fidelity reconstruction ability was verified through comparison with analytical solutions. To validate the applicability of the framework under multi-material conditions, copper, aluminum, steel, and titanium were selected as representative metals for analysis. The temperature distribution characteristics of the pin fin were compared under different thermal conductivities. The results indicate a significant positive correlation between the pin-fin tip temperature and the material thermal conductivity, further proving the framework’s ability to accurately capture material-dependent thermal effects. The main contributions of this paper include: (i) the proposal of a generalizable strategy for transforming boundary value problems and solving initial value problems; (ii) a systematic evaluation of the numerical accuracy and stability of pin-fin thermal modeling in multi-material systems; (iii) providing a quantitative engineering basis for multi-material pin-fin heat dissipation designs under high heat flux density conditions. Overall, this study not only deepens the understanding of the heat dissipation mechanisms of pin fins under different material conditions but also provides a feasible algorithmic approach for high-precision numerical modeling of microscale thermal systems.

2.Related work

2.1 Development and Application of Classical Numerical Solvers
Early heat conduction/heat transfer numerical solvers often employed basic integration methods such as explicit Euler, implicit trapezoidal, and classical fourth-order Runge-Kutta methods, which were widely used in engineering due to their simplicity and low cost. However, under high heat flux density and mixed boundary conditions, low-order methods often accumulate truncation errors rapidly unless the step size is significantly reduced, thereby sacrificing efficiency. Modern numerical analysis also indicates that commonly used implicit formats (such as backward Euler and trapezoidal methods) can be equivalently represented by Runge-Kutta structures, which helps in understanding their stability characteristics in (semi)stiff scenarios [3]. In practical applications, higher-order explicit methods with embedded error estimation and adaptive step size, such as the Dormand-Prince family and its implementation in MATLAB ode45, have become the “default tool” for solving non-stiff ordinary differential equations. However, without proper boundary value problem reconstruction, their direct application may still be inefficient [4][5].

2.2 Numerical Strategies for Boundary Value Problems
Two-point boundary value problems are challenging to solve due to the distribution of boundary conditions at different locations. The shooting method transforms BVPs into a series of IVPs, iteratively adjusting the unknown initial derivatives to satisfy the terminal boundary. In terms of root-finding update strategies, the secant method is particularly practical because it does not require derivative information and, under well-behaved mappings, has a convergence rate similar to that of Newton’s method. This makes it especially suitable for scenarios where sensitivity is difficult to obtain or costly [6][7]. Comparative studies have also placed the shooting method alongside finite differences, collocation/finite element methods, highlighting the trade-offs between robustness and implementation complexity [8]. However, in the field of electronic heat dissipation, there is still a lack of systematic validation for the shooting + secant method applied to micro-scale, multi-material pin fins with mixed boundary conditions, which has motivated the proposal of the unified framework in this paper.

2.3 Advantages of High-Order Explicit Integration Methods
In scenarios demanding high accuracy, embedded high-order explicit methods can achieve excellent error control at relatively low cost. The Dormand-Prince (RK5(4)) family is well known for its embedded pair and FSAL (First Same As Last) strategy, requiring only six function evaluations to perform adaptive step-size adjustment [9][10]. Its widespread use in solvers such as MATLAB ode45 is attributed to reliable local truncation error estimation, dense output interpolation, and robust performance across a wide range of smooth problems [4][11]. However, in boundary value problems (BVPs) such as pin-fin conduction—where geometric scales, material contrasts, and boundary types are strongly coupled—the synergistic effect between high-order explicit integration and slope updating (e.g., secant-based shooting) remains insufficiently characterized, particularly under multi-material parameter sweeps. Establishing such coupling is crucial for simultaneously achieving accuracy and efficiency in design-oriented studies.

2.4 Current Status of Multi-Material Heat Dissipation Modeling
Research on pin-fin cooling has long focused on geometric optimization and flow-heat transfer coupling (e.g., array spacing, inline vs. staggered arrangements, and pressure drop trade-offs), with numerous experimental and numerical studies reporting heat transfer enhancement mechanisms [12][13]. Analytical pin-fin models remain central to conduction modeling under convective or mixed boundary conditions, and in recent years have been extended to include radiation and temperature-dependent properties [14][15]. However, systematic algorithm-level comparisons of solver accuracy and stability across multiple materials (e.g., copper, aluminum, steel, and titanium) remain relatively scarce; many studies tend to fix the material while varying only geometry or flow conditions. Consequently, it is necessary to evaluate a unified numerical framework that can maintain both accuracy and robustness under conditions of large thermal conductivity contrasts and multiple boundary constraints.

2.5 Differentiated Contributions of This Study
To address the aforementioned gaps, this study integrates the “secant-assisted shooting method” with a “low-order Euler predictor + high-order Dormand-Prince corrector,” forming a unified BVPs→IVPs workflow tailored for pin-fin problems. The framework combines rapid slope identification without derivatives with high-fidelity temperature field reconstruction and is stress-tested across multiple representative heat dissipation materials. Unlike prior studies that primarily emphasized geometry or turbulence enhancement, this work highlights the solver’s transferability under varying materials and mixed boundary conditions, aligning with the recent need for stable numerical kernels to enable fine-grained optimization of perforated/microstructured pin-fins [16][4]. This establishes a closed loop linking materials, numerics, and design, thereby providing a reusable toolchain for multi-material pin-fin analysis.

3.Method

3.1 Physical and Mathematical Modeling
This study considers the steady-state one-dimensional heat conduction process in pin fins, assuming that the material is isotropic, with constant thermophysical properties (such as thermal conductivity and specific heat capacity), and that the convective heat transfer coefficient is uniformly distributed over the fin surface. Based on the principle of energy conservation and Fourier’s law, the heat transfer process can be described by the following second-order ordinary differential equation:

1) Steady state heat conduction equation:

(1)

where is the characteristic parameter of the equation, is the perimeter of the pin fin,is the cross-sectional area of the fin,is the thermal conductivity of the material,is the convective heat transfer coefficient, andis the ambient temperature.

2) boundary conditions:

The temperature at the base is known, denoted . The end of the pin fin is assumed to be at a boundary condition, where the temperature gradient is zero, i.e. . Here, is the length of the pin fin.

3.2 Initial Value Problem Reformulation

To transform the second-order boundary value problem into a first-order system, we introduce new state variables resulting in the following set of first-order equations:

(2)

where denotes the temperature and denotes its derivative. By choosing an initial slope and determining it through the shooting method (see §3.3), the boundary conditions can be satisfied.

3.3 Shooting Method and Secant Update

1) Secant Iteration: The shooting method transforms the boundary value problem into an initial value problem, with the unknown initial derivative treated as an iterative variable. Numerical integration is performed from the known boundaryto the adiabatic end, yielding a residual function for the terminal temperature gradient. The value of the residual function, representing the deviation from the boundary condition, is then used to update the guess for the initial derivative. By employing the secant method to construct the update formula, the complexity of derivative evaluation is avoided, thereby improving computational efficiency and reformulating the complex boundary value problem into an initial value problem.

(3)

Wheredenotes the residual function, the calculation formula is, defined as the deviation between the numerical solution at the fin tip and the adiabatic boundary condition.

• Convergence Criterion:The convergence condition is satisfied when the difference between two successive estimates of the initial derivative is smaller than a prescribed tolerance, typically , i.e.:

(4)

Whereis the iterative variable,. In this study, the convergence tolerance for the initial derivative in the shooting-secant iteration is set to,  since this threshold provides a practical compromise between numerical analysis and engineering requirements. On the one hand, the machine precision of double-precision floating-point arithmetic is on the order of ; further tightening the tolerance would significantly increase the iteration count without yielding perceivable improvements in physical accuracy, while amplifying the impact of round-off errors [17]. On the other hand, the default absolute error tolerance of mainstream ODE solvers is also set at, which typically corresponds to terminal errors below the milli-Kelvin level in temperature fields, sufficient to satisfy engineering criteria and robustness requirements [18][4]. Considering the propagation of truncation and round-off errors, the superlinear convergence property of the secant method, and alignment with the error control strategies of subsequent numerical integrators (e.g., Dormand-Prince/ode45), the tolerance ofensures both solution accuracy and compliance with boundary conditions while avoiding unnecessary computational overhead [19][20].

3.4 Numerical Integrator Design

To achieve numerical solutions, two types of numerical integrators are introduced in this study: the low-order explicit Euler method and higher-order Runge-Kutta (RK5) methods.

• Low-Order Baseline: The explicit Euler method is employed as a low-order integration benchmark to verify stability under a simple solver framework.

(5)

where the step size , ensuring comparability of results across different integration methods.

2) High-Order Main Scheme: A fifth-order Runge-Kutta method based on Dormand-Prince coefficients is implemented to obtain high-accuracy temperature distributions. At each step, the high-order integrator requires multiple evaluations of the conduction equation’s slope information, thereby improving global accuracy and mitigating the accumulation of truncation errors.

(6)   

(7)

(8)

(9)

(10)

(11)

(12)

) (13)

(14)

(15)

(16)

(17)

(18)

(19)

Formulas (6-12) and (13-19) reflect the core idea of the fifth-order Runge-Kutta method for solving the first-order system in the interval multiple slopes …,  , are constructed and weighted to obtain a higher precision of .

3.5 Analytical Solution for Consistency Check
To verify the correctness of the numerical method, an analytical solution is obtained using the symbolic computation tool Matlab’s dsolve, which serves as a benchmark for the numerical solution. For constant material properties and simple geometric shapes, the general form of the analytical solution is:

(20)

By applying the boundary conditionsand , the constants are solved respectively:

(21)

Here,is the temperature of the fin base,is the ambient temperature,, and 𝑚 is the characteristic parameter.

4. Experimental Design and Results Analysis

4.1 Objectives and Overall Approach
The experimental component of this study aims to validate the applicability, accuracy, and stability of the proposed numerical framework-combining the shooting method, secant iteration, and high-order explicit integration—for multi-material pin-fin heat dissipation modeling. The focus is placed on comparing the temperature distribution characteristics of representative metallic materials with different thermal conductivities under identical geometric and boundary conditions, and on evaluating the consistency and robustness of the proposed method across a wide material parameter space.

4.2 Materials and Thermophysical Parameter Settings
Four common metals—copper (Cu), aluminum (Al), steel (Steel), and titanium (Ti)—are selected as representative materials, covering a typical range from high to low thermal conductivity. For each material, thermal conductivity, density, specific heat capacity, convective heat transfer coefficient, and ambient temperature are specified according to standard engineering reference data and summarized in Table 1. All environmental variables are kept identical across different materials to ensure representativeness of the simulation conditions. A uniform cross-sectional pin-fin model is employed to eliminate geometric influences, with identical fin length, cross-sectional area, and perimeter assigned to all materials. The primary thermophysical parameters and environmental conditions for the four metals are listed below.

Table1. Main physical and chemical properties of four metal materials and environmental conditions

4.3 Numerical Solver Configuration

1) Integration step size: fixed at to ensure comparability of results across different methods;

2) Convergence tolerance: set toto guarantee stable convergence of the shooting-secant iteration;

3) Numerical methods: both the explicit Euler method and the fifth-order Runge-Kutta method (Dormand-Prince) are employed for numerical integration, with comparisons made against the analytical solution to assess accuracy and stability.

4.4 Main Results and Comparative Analysis

1) In the simulations, high-thermal-conductivity materials (e.g., copper and aluminum) exhibit gradual temperature decay curves, indicating that heat is more effectively conducted within the material. As a result, the tip temperatures are significantly higher than those of low-thermal-conductivity materials (e.g., steel and titanium). In particular, due to their relatively large thermal conductivities, copper and aluminum display slower temperature variation, allowing heat to distribute more uniformly over longer distances. In contrast, steel and titanium, as low-conductivity materials, show steep temperature drops along the fin length, leading to tip temperatures much lower than those of the former. This phenomenon is directly linked to the thermal conductivity of the material and is consistent with physical expectations—that materials with higher thermal conductivity can transfer heat more effectively and thus maintain higher tip temperatures.

Fig.1. Temperature distribution curve of copper, aluminum, steel and titanium materials in the length direction under steady state heat conduction

As shown in Fig.1, the temperature distribution curves of different materials exhibit pronounced differences. Copper and aluminum display relatively gradual temperature decay curves, maintaining higher temperatures at the fin tip, whereas steel and titanium show rapid temperature drops over shorter conduction distances. This result clearly demonstrates the dominant role of thermal conductivity in steady-state heat transfer processes.

• By comparing the tip temperatures of the four materials (at), a consistent ranking is obtained: copper > aluminum > steel > titanium. The temperature differences are positively correlated with the thermal conductivity, indicating that materials with higher conductivity are able to transfer heat over longer distances, thereby resulting in higher tip temperatures. Copper, with its superior thermal conductivity, achieves the highest tip temperature, while titanium, with its relatively low conductivity, shows the lowest. This ranking validates the significant influence of thermal conductivity on the temperature distribution.

Fig.2. Comparison of temperature of copper, aluminum, steel and titanium at the end position (x=0.04).

Fig. 2 compares the tip temperatures of the four materials. The results show that copper and aluminum maintain relatively high temperatures at the fin tip, whereas steel and titanium exhibit significantly lower values. This outcome is not only consistent with the relative magnitudes of their theoretical thermal conductivities but also further validates the reliability of the numerical model in predicting material heat transfer performance.

• To verify the accuracy of the numerical solution methods, the results of the fifth-order Runge-Kutta method were compared against the analytical solution. The findings indicate that the fifth-order Runge-Kutta method matches the analytical solution with excellent agreement across the entire domain, with negligible errors, demonstrating its high accuracy and stability in handling heat conduction problems. In contrast, while the Euler method captures the overall trend, it shows slight deviations in the tip region (x→L), which can be attributed to approximation errors inherent in the method when using relatively large step sizes.

To further quantify the accuracy differences between the two numerical methods, the root-mean-square error (RMSE) relative to the analytical solution was calculated. By comparing the RMSE curves, the error levels of different numerical methods under various material conditions were assessed.

(22)

Where N is the total number of sampling points,is the temperature value of the numerical method (Euler method or Runge-Kutta method) at the position. is the temperature value of the analytical solution at the same position.

Fig.3. Comparison of root mean square error (RMSE) between numerical method and analytical solution at along (x=0→0.04)

Table2. Two numerical solution RMSE errors of four materials at the end of the pin-fin

As shown in Fig.3 and Table 2, the fifth-order Runge-Kutta method exhibits negligible deviations throughout the entire heat conduction process, whereas the explicit Euler method shows noticeable errors near the fin tip, which become particularly pronounced under high-thermal-conductivity material conditions.

5. Summary and Outlook

5.1 Framework and Innovations

This study proposes and validates an efficient numerical framework for solving steady-state heat conduction boundary value problems. By reformulating the original boundary value problem into an equivalent initial value problem, the method integrates the secant approach for rapid iteration of the initial derivative and employs both the explicit Euler method and the high-order Dormand-Prince fifth-order Runge-Kutta method for integration, achieving a balance between computational efficiency and numerical accuracy. In parallel, an analytical solution was established as a benchmark, providing a rigorous validation mechanism for the numerical results. Accordingly, this study directly addresses the critical challenge highlighted in the introduction: how to achieve reliable, scalable, and efficient solutions for boundary-driven heat conduction problems in engineering thermal design.

5.2 Numerical Performance and Applicability

Numerical experiments conducted under four representative materials—copper, aluminum, steel, and titanium—demonstrated that the proposed method exhibits stable and consistent convergence across a wide range of thermal conductivities. The results show that in all computational cases the number of convergence iterations remained below ten, underscoring the efficiency of the framework. Comparative analyses further confirmed that the high-order Runge-Kutta method significantly reduces global errors under identical step sizes, while the Euler method, despite its limited accuracy, still serves as a fast and practical tool for preliminary computations. Through systematic comparisons with the analytical solution, this study ensures that the results are not only computationally efficient but also physically consistent and reliable.

5.3 Physical Insights and Engineering Implications

Beyond methodological contributions, the simulation results of this study reveal fundamental physical insights into heat conduction. In particular, the tip temperature is shown to exhibit a strong positive correlation with the material’s thermal conductivity, explaining why high-conductivity metals (e.g., copper and aluminum) can retain elevated temperatures at distal boundaries. This finding holds both academic significance and practical value, providing quantitative guidance for material selection, geometric optimization, and performance prediction in engineering heat sinks, especially pin-fin structures. Furthermore, the proposed framework demonstrates strong generality, suggesting its potential applicability to other boundary-driven heat conduction problems, thereby bridging numerical mathematics and engineering thermophysics.

5.4 Limitations and Future Work

Nevertheless, this study has certain limitations that warrant further exploration. First, the integration process adopted a fixed step size, which restricts the dynamic balance between accuracy and efficiency. Future work could incorporate adaptive step-size control strategies, particularly those based on embedded high-order methods, to achieve superior computational performance. Second, although the secant method exhibited high efficiency in this study, its convergence rate may deteriorate in strongly nonlinear systems or in cases with drastic parameter variations. The integration of modern numerical tools such as automatic differentiation may further accelerate convergence. Finally, the applicability of this framework has so far been restricted to one-dimensional steady-state heat conduction. Extending the approach to transient heat conduction, multi-dimensional geometries, and multi-physics coupling problems represents a natural direction for future research, especially given the increasing complexity of modern engineering thermal management systems.

5.5 Overall Conclusion

In summary, the numerical framework proposed in this study combines theoretical rigor with engineering practicality, offering a powerful and extensible tool for thermal modeling at micro- and mesoscopic scales. Through validation against analytical solutions and demonstrations under multi-material conditions, the method has proven its stability, accuracy, and efficiency in addressing boundary-driven heat conduction problems. Looking forward, as the current limitations are progressively overcome and the framework is extended to complex geometries and coupled-field problems, this line of research is expected to play an increasingly profound role in advancing both the scientific understanding of heat conduction and the practical design of engineering thermal systems.

Reference

[1] Wang Z, Dong R, Ye R, Singh S S K, Wu S, Chen C. A review of thermal performance of 3D stacked chips [J]. International Journal of Heat and Mass Transfer, 2024, 235: 126212.

[2] Garimella S V, Fleischer A S, Murthy J Y, Keshavarzi A, Prasher R, Patel C, et al. Thermal challenges in next-generation electronic systems [J]. IEEE Transactions on Components and Packaging Technologies, 2008, 31(4): 801-815.

[3] Atkinson K, Han W. Theoretical Numerical Analysis: A Functional Analysis Framework [M]. 3rd ed. New York, NY, USA: Springer, 2009.

[4] MathWorks. ode45 — Solve nonstiff differential equations [OL]. Documentation, 2024.

[5] Dormand J R, Prince P J. A family of embedded Runge-Kutta formulae [J]. Journal of Computational and Applied Mathematics, 1980, 6(1): 19-26.

[6] Iron P. Numerical solutions of boundary value problems: Shooting—Secant method [M]. Lecture Notes, Dalhousie University, 2012.

[7] Brigham Young University. The shooting method for boundary value problems [OL]. Notes, 2017.

[8] Heckbert P. Boundary value problems [M]. In: 15-859B: Introduction to Scientific Computing. Carnegie Mellon University, 2000.

[9] Dormand J R, Prince P J. New Runge-Kutta algorithms for numerical simulation in dynamical astronomy [J]. Celestial Mechanics, 1978, 18: 223-232.

[10] Moler C B. Ordinary differential equation solvers: ODE23 and ODE45 [EB]. Cleve’s Corner, MathWorks Blog, 2014.

[11] Calvo M, Montijano J I, Rández L. A fifth-order interpolant for the Dormand and Prince Runge-Kutta method [J]. Journal of Computational and Applied Mathematics, 1990, 29(1): 91-100.

[12] Lawson S A, Beér O, Hamdan M O, Arts T. Heat transfer from multiple row arrays of low aspect ratio pin fins [J]. International Journal of Heat and Mass Transfer, 2011, 54: 4099-4109.

[13] Comparison of pin-surface heat transfer in arrays of oblong pin fins [J]. Journal of Turbomachinery, ASME, 2013.

[14] MIT Unified Engineering Notes. 18.2 Heat transfer from a fin [OL]. Massachusetts Institute of Technology, 2010s.

[15] Alhejaili W, Almurayshid H. Analytical solution for temperature equation of a fin including convection and radiation with internal heat generation [J]. Chemical Physics Letters, 2022, 790: 139404.

[16] Ehsani H, Akbari M A, Nazari M. Investigating thermal performance enhancement in perforated pin-fin heat sinks [J]. Case Studies in Thermal Engineering, 2024, 50: 104029.

[17] Higham N J. Accuracy and Stability of Numerical Algorithms [M]. 2nd ed. Philadelphia, PA, USA: SIAM, 2002.

[18] Burden R L, Faires J D. Numerical Analysis [M]. 9th ed. Boston, MA, USA: Brooks/Cole, 2010.

[19] MathWorks Community. Default tolerance of ode45 [EB].

[20] Wikipedia. Machine epsilon [OL].