# Generalized Laplace-Type Transform Method for Solving Multilayer Diffusion Problems

Authors: Mohamed Akel; Hillal M. Elshehabey; Ragaa Ahmed
Journal: Journal of Function Spaces
 (2022)

---

## Abstract

Multilayer diffusion problems have found significant importance that they arise in many medical, environmental, and industrial applications of heat and mass transfer. In this article, we study the solvability of a one-dimensional nonhomogeneous multilayer diffusion problem. A new generalized Laplace-type integral transform is used, namely, theMρ,m-transform. First, we reduce the nonhomogeneous multilayer diffusion problem into a sequence of one-layer diffusion problems including time-varying given functions, followed by solving a general nonhomogeneous one-layer diffusion problem via the Mρ,m-transform. Hence, by means of general interface conditions, a renewal equations’ system is determined. Finally, the Mρ,m-transform and its analytic inverse are used to obtain an explicit solution to the renewal equations’ system. Our results are of general attractiveness and comprise a number of previous works as special cases.


---

## Body

## 1. Introduction

The multilayer diffusion problems are typical models for a variety of solute transport phenomena in layered permeable media, such as advection, dispersion, and reaction diffusions [1–10]. These problems have had their importance due to their natural prevalence in a remarkably large number of applications such as chamber-based gas flux measurements [11], contamination and decontamination in permeable media [6, 12], drug-eluting stent [13, 14], drug absorption [15, 16], moisture propagation in woven fabric composites [17], permeability of the skin [18], and wool-washing [19]. Further applications have been considered in [20, 21].As epidemiological models, reaction-diffusion problems are widely used to model and analyze the spread of diseases such as the global COVID-19 pandemic caused by SARS-CoV-2. These models describe the spatiotemporal prevalence of the viral pandemic and apprehend the dynamics depending on human habits and geographical features. The models estimate a qualitative harmony between the simulated prediction of the local spatiotemporal spread of a pandemic and the epidemiological collected datum (see [22, 23]). These data-driven emulations can essentially inform the respective authorities to purpose efficient pandemic-arresting measures and foresee the geographical distribution of vital medical resources. Moreover, such studies explore alternate scenarios for the repose of lockdown restrictions based on the local inhabitance densities and the qualitative dynamics of the infection. For more applications, one can refer, e.g., to [24, 25].Although the numerical methods are usually applied to solve the diffusion problems, especially in the heterogeneous permeable media, the analytic solutions, when available, are characterized by their exactness and continuity in space and time. In the context of obtaining numerical solutions for such models, we refer to the following references [26–30]. In this work, we focus on analytic solutions of certain nonhomogeneous diffusion problems in multilayer permeable media. Here, the retardation factors are assumed to be constant, the dispersion coefficients vary across layers, but being constants within each layer, and the free terms are (arbitrary) time-varying functions.Analytic and semianalytic solutions of multilayer diffusion problems are developed by using the integral transforms [6, 25, 31–40]. Applying Laplace transforms, to solve multilayer diffusion problems, has advantages as an applicable tool in handling different types of boundary conditions and averts solving complicated transcendental equations as demanded by eigenfunction expansion methods. Further works involving the Laplace transform have studied the permeable layered reaction diffusion problem in [41, 42]. Solutions obtained in these works are restricted to two layers as well as obtaining the inverse Laplace transform numerically. In the same context, generalized integral transform techniques, for short GITT, are well-established hybrid approaches for solving diffusion and convection-diffusion problems, in which hybrid refers to the combination of classical analytical methods with modern computational tools aimed at accurate, robust, and low-cost solutions [43–47]. In the current work, we aim to extend, generalize, and merge results in [31, 33, 38, 40, 42] to solve certain nonhomogeneous diffusion problems in one-dimensional n-layered media. We use a new generalized integral transform recently introduced in [48]. The obtained solutions are applicable to more general linear nonhomogeneous diffusion equations, finite media consisting of arbitrary many layers, continuity and dispersive flow at the contact interfaces between sequal layers, and transitory boundary conditions of the arbitrary type at the inlet and outlet. To the best knowledge of the authors, analytical solutions verifying all the above mentioned conditions have not been previously reported in the literature which strongly motivates this current work.In the remaining part of this introductory section, in Subsection1.1, the multilayer diffusion problem is described, and then, it is reformulated as a sequence of one-layer diffusion problems having boundary conditions including given time-depending functions. Basic properties for theMρ,m-transform that will be needed in this work are stated in Subsection 1.2. The remaining sections are constructed as follows: in Section 2, we discuss the solvability of a general linear nonhomogeneous one-layer diffusion problem with arbitrary time-varying data, using the Mρ,m-transform. Section 3 is devoted to our main multilayer diffusion problem, where in Subsection 3.1, we solve a two-layer problem to shed light on the basic idea by considering this simple case. Further, in Subsection 3.2, we return to benefit from the results obtained in Section 2 and Subsection 3.1 to solve the main multilayer diffusion problems (2)–(8) (see Subsection 1.1 below).


### 1.1. Mathematical Modeling for Nonhomogeneousn-Layer Diffusion Systems

A one-dimensional diffusion problem in ann-layered permeable medium is set out as follows. Let
(1)α=x0<x1<⋯<xn−1<xn=β,be a finite partition of the interval α,β. In each subinterval xj−1,xj, with j=1,2,⋯,n, the component function φjx,t satisfies the partial differential equation (PDE)
(2)∂φj∂t=dj∂2φj∂x2+λt,τrjx,t,x∈xj−1,xj,t,τ>0,where dj≥0, for all 1≤j≤n, are the diffusion coefficients and
(3)λt,τ=tmτm+τm−ρ,t≥0,τ>0,with m∈ℤ+=1,2,3,⋯,ρ∈ℂ,Reρ>0. Here, the function-term λt,τrjx,t physically means the external source term that could be applied to the diffusion equation with rjx,t depends on time and space while the other factor of the source term, i.e., λ, depends only on time. This last factor could be, for instance, a periodic-time magnetic source.The initial conditions (ICs) are assumed as(4)φjx,0=ηjx,x∈xj−1,xj,1≤j≤n.The boundary conditions (BCs) are posited as(i)
The outer BCs (at the inletx=α and the outlet x=β) are general Robin boundary conditions as(5)ıφ1α,t+ι∂φ1∂xα,t=λt,τζt,(6)ℓφnβ,t+l∂φn∂xβ,t=λt,τξt,for allt≥0, with ı,ι,ℓ, and l are constants satisfying ∣ı∣+∣ι∣>0,∣ℓ∣+∣l∣>0.
(ii)
The inner BCs (the interface conditions) are(7)φjxj,t=Λjφj+1xj,t,(8)νjφjxj,t+μj∂φj∂xxj,t=νj+1φj+1xj,t+μj+1∂φj+1∂xxj,t,for all t≥0, with ∣νj∣+∣μj∣>0 for all j=1,2,⋯,n−1.For appropriate given functionsη1,⋯,ηn,ζ, and ξ, we are going to find an analytic solution of the problems (2)–(8) using the Mρ,m-generalized integral transform, introduced recently in [48]. Problems (2)–(8) can be reduced into the following sequence of one-layer diffusion problems.
(i)
In the inlet layer, i.e.,x∈x0,x1,(9)∂φ1∂t=d1∂2φ1∂x2+λt,τr1x,t,x∈x0,x1,t,τ>0,φ1x,0=η1x,x∈x0,x1,ιφ1x0,t+ι∂φ1∂xx0,t=λt,τζ1t,t≥0,τ>0,ν1φ1x1,t+μ1∂φ1∂xx1,t=λt,τξ1t,t≥0,τ>0.(ii)
In the interior layers, i.e.,x∈xj−1,xj,2≤j≤n−1,(10)∂φj∂t=dj∂2φj∂x2+λt,τrjx,t,x∈xj−1,xj,t,τ>0,φjx,0=ηjx,x∈xj−1,xj,νjφjxj−1,t+μj∂φj∂xxj−1,t=λt,τζjt,t≥0,τ>0,νjφjxj,t+μj∂φj∂xxj,t=λt,τξjt,t≥0,τ>0.(iii)
In the outlet layer, i.e.,x∈xn−1,xn,(11)∂φn∂t=dn∂2φn∂x2+λt,τrnx,t,x∈xn−1,xn,t,τ>0,φnx,0=ηnx,x∈xn−1,xn,νnφnxn−1,t+μn∂φn∂xxn−1,t=λt,τζnt,t≥0,τ>0,ℓφnxn,t+l∂φn∂xxn,t=λt,τξnt,t≥0,τ>0.Remark 1.
Each of the initial boundary value problems (9)–(11) is a case of the one-layer nonhomogeneous diffusion problem that will be discussed in Section 2 below.Now, in view of the inner boundary conditions (7) and (8), the time-varying functions ζj and ξj for all 2≤j≤n are subject to
(12)ζjt=ξj−1t,2≤j≤n,so that
(13)ζjt=ζt,j=1,ξj−1t,2≤j≤n,ξjt=νjφjxj,t+μj∂φj∂xxj,tλ−1t,τ,1≤j≤n−1,ξt,j=n.While the outer boundary dataζ1t=ζt and ξnt=ξt are given in (5) and (6), respectively, the functions ζj2≤j≤n can be determined once we specify the functions ξj1≤j≤n−1. Hence, we have to find ξj,1≤j≤n−1. To do so, we should use the first matching condition (7).


### 1.2. Srivastava-Luo-Raina Generalized Integral Transform

In [48], Srivastava et al. introduced the following generalized integral transform:
(14)Mρ,mφts,τ=∫0∞e−stφτttm+τmρdt,for a continuous (or piecewise continuous) function φ on 0,∞, where ρ∈ℂ;Reρ≥0;m∈ℤ+,s>0 is the transform variable and τ>0 is a parameter. The basic properties of the Mρ,m-transform are given in [48]. Next, we recall some of these properties, which are needed in the present work. Indeed, as introduced in [48] the Mρ,m-transform is closely related with the well-known integral transforms, the Laplace, natural, and Sumudu transforms. The Laplace transform is defined by
(15)Lφts=∫0∞e−stφtdt,Res>0.So, from (14) and (15), we have the following duality relations:
(16)Lφts=M0,mφts,1,Res>0,Mρ,mφts,τ=Lφτttm+τmρs,s,τ>0,Mρ,mφts,τ=1τLφttm/τm+τmρsτ,s,τ>0,(17)Mρ,mtmτm+τmρφts,τ=Lφτts,s,τ>0.Settingρ=0 in (14), we recover the natural transform defined as (see [49, 50])
(18)ℕφts,τ=∫0∞e−stφτtdt,s>0,τ>0.Thus, we have the followingMρ,m‐ℕ-transform duality
(19)ℕφts,τ=M0,mφts,τ,(20)Mρ,mφts,τ=ℕφttm/τm+τmρs,τ,s>0,τ>0,(21)Mρ,mtmτm+τmρφts,τ=ℕφts,τ;s>0,τ>0.The Sumudu transform is defined by [51–53]
(22)Sφtτ=∫0∞e−tφτtdt,τ>0.Thus,(23)Sφtτ=M0,mφt0,τ,τ>0,(24)Mρ,mφts,τ=1sSφttm/τm+τmρτs,s,τ>0.Based on these dualities of theMρ,m-transform (14) and these well-known integral transforms, it seems to be interesting to apply the Mρ,m-transform (14) in solving a variety of boundary and initial-boundary value problems. In this context, we recall the following results [48]:
(i)
Letφnt be the nth-order t-derivative of the function φt and φt≤Ket/γ with K>0, γ>0. Then,(25)Mρ,mtmτm+τmρφnts,τ=snτnℕφts,τ−∑k=0n−1skτk+1φn−k−10,where ℕφts,τ is defined by (18). Using the duality (21) in (25), we find
(26)ℕφnts,τ=snτnℕφts,τ−∑k=0n−1skτk+1φn−k−10,n=0,1,⋯.(ii)
Again, using the dualities stated before a convolution formula for theMρ,m-transform (14) can be obtained as follows. Here, the convolution for the Laplace transform will be considered; that is, for the functions φ and ψ, the convolution formula is given as(27)φ∗ψt=∫0tφxψt−xdx=∫0tφt−xψxdx.IfΦs,τ=Mρ,mφts,τ and Ψs,τ=Mρ,mψts,τ, thenτΦs,τΨs,τ=τ∫0∞e−st1φτt1/t1m+τmρdt1∫0∞e−st2ψτt2/t2m+τmρdt2=τ∫0∞∫0∞e−st1+t2φ~τt1ψ~τt2dt1dt2,where
(28)φ~t=φttm/τm+τmρ,ψ~t=ψttm/τm+τmρ.Settingt1+t2=t in the last equality, one gets
(29)τΦs,τΨs,τ=τ∫0∞∫t2∞e−stφ~τt−t2ψ~τt2dtdt2=τ∫0∞e−stdt∫0tφ~τt−t2ψ~τt2dt2.Here, changing of the integral order is used. Thus, using the duality of the Mρ,m and ℕ transforms (see (20)), we find
(30)τΦs,τΨs,τ=ℕφ~∗ψ~s,τ.Remark 2.
If we putρ=0 in (30), the case being interesting later in our work, then we get
(31)τℕφts,τℕψts,τ=ℕφ∗ψts,τ.(iii)
Once again, ...again, using the dualities stated before an inversion formula of theMρ,m-transform (14) is given (see Theorem 4.1 of [48]) as(32)φt=tmτm+τmρL−1Mρ,mφts,ttτ=12πitmτm+τmρ∫c−i∞c+i∞est/τMρ,mφts,τds,c,τ>0,as long as the integral converges absolutely. In case, when ρ=0 one obtains the following inversion formula of the natural transform (see Theorem 5.3 of [49])
(33)φt=12πi∫c−i∞c+i∞est/τℕφts,τds,c,τ>0.The residue theorem (see, e.g., [54]) is usually used to calculate the contour integrals in (32) and (33).


## 1.1. Mathematical Modeling for Nonhomogeneousn-Layer Diffusion Systems

A one-dimensional diffusion problem in ann-layered permeable medium is set out as follows. Let
(1)α=x0<x1<⋯<xn−1<xn=β,be a finite partition of the interval α,β. In each subinterval xj−1,xj, with j=1,2,⋯,n, the component function φjx,t satisfies the partial differential equation (PDE)
(2)∂φj∂t=dj∂2φj∂x2+λt,τrjx,t,x∈xj−1,xj,t,τ>0,where dj≥0, for all 1≤j≤n, are the diffusion coefficients and
(3)λt,τ=tmτm+τm−ρ,t≥0,τ>0,with m∈ℤ+=1,2,3,⋯,ρ∈ℂ,Reρ>0. Here, the function-term λt,τrjx,t physically means the external source term that could be applied to the diffusion equation with rjx,t depends on time and space while the other factor of the source term, i.e., λ, depends only on time. This last factor could be, for instance, a periodic-time magnetic source.The initial conditions (ICs) are assumed as(4)φjx,0=ηjx,x∈xj−1,xj,1≤j≤n.The boundary conditions (BCs) are posited as(i)
The outer BCs (at the inletx=α and the outlet x=β) are general Robin boundary conditions as(5)ıφ1α,t+ι∂φ1∂xα,t=λt,τζt,(6)ℓφnβ,t+l∂φn∂xβ,t=λt,τξt,for allt≥0, with ı,ι,ℓ, and l are constants satisfying ∣ı∣+∣ι∣>0,∣ℓ∣+∣l∣>0.
(ii)
The inner BCs (the interface conditions) are(7)φjxj,t=Λjφj+1xj,t,(8)νjφjxj,t+μj∂φj∂xxj,t=νj+1φj+1xj,t+μj+1∂φj+1∂xxj,t,for all t≥0, with ∣νj∣+∣μj∣>0 for all j=1,2,⋯,n−1.For appropriate given functionsη1,⋯,ηn,ζ, and ξ, we are going to find an analytic solution of the problems (2)–(8) using the Mρ,m-generalized integral transform, introduced recently in [48]. Problems (2)–(8) can be reduced into the following sequence of one-layer diffusion problems.
(i)
In the inlet layer, i.e.,x∈x0,x1,(9)∂φ1∂t=d1∂2φ1∂x2+λt,τr1x,t,x∈x0,x1,t,τ>0,φ1x,0=η1x,x∈x0,x1,ιφ1x0,t+ι∂φ1∂xx0,t=λt,τζ1t,t≥0,τ>0,ν1φ1x1,t+μ1∂φ1∂xx1,t=λt,τξ1t,t≥0,τ>0.(ii)
In the interior layers, i.e.,x∈xj−1,xj,2≤j≤n−1,(10)∂φj∂t=dj∂2φj∂x2+λt,τrjx,t,x∈xj−1,xj,t,τ>0,φjx,0=ηjx,x∈xj−1,xj,νjφjxj−1,t+μj∂φj∂xxj−1,t=λt,τζjt,t≥0,τ>0,νjφjxj,t+μj∂φj∂xxj,t=λt,τξjt,t≥0,τ>0.(iii)
In the outlet layer, i.e.,x∈xn−1,xn,(11)∂φn∂t=dn∂2φn∂x2+λt,τrnx,t,x∈xn−1,xn,t,τ>0,φnx,0=ηnx,x∈xn−1,xn,νnφnxn−1,t+μn∂φn∂xxn−1,t=λt,τζnt,t≥0,τ>0,ℓφnxn,t+l∂φn∂xxn,t=λt,τξnt,t≥0,τ>0.Remark 1.
Each of the initial boundary value problems (9)–(11) is a case of the one-layer nonhomogeneous diffusion problem that will be discussed in Section 2 below.Now, in view of the inner boundary conditions (7) and (8), the time-varying functions ζj and ξj for all 2≤j≤n are subject to
(12)ζjt=ξj−1t,2≤j≤n,so that
(13)ζjt=ζt,j=1,ξj−1t,2≤j≤n,ξjt=νjφjxj,t+μj∂φj∂xxj,tλ−1t,τ,1≤j≤n−1,ξt,j=n.While the outer boundary dataζ1t=ζt and ξnt=ξt are given in (5) and (6), respectively, the functions ζj2≤j≤n can be determined once we specify the functions ξj1≤j≤n−1. Hence, we have to find ξj,1≤j≤n−1. To do so, we should use the first matching condition (7).


## 1.2. Srivastava-Luo-Raina Generalized Integral Transform

In [48], Srivastava et al. introduced the following generalized integral transform:
(14)Mρ,mφts,τ=∫0∞e−stφτttm+τmρdt,for a continuous (or piecewise continuous) function φ on 0,∞, where ρ∈ℂ;Reρ≥0;m∈ℤ+,s>0 is the transform variable and τ>0 is a parameter. The basic properties of the Mρ,m-transform are given in [48]. Next, we recall some of these properties, which are needed in the present work. Indeed, as introduced in [48] the Mρ,m-transform is closely related with the well-known integral transforms, the Laplace, natural, and Sumudu transforms. The Laplace transform is defined by
(15)Lφts=∫0∞e−stφtdt,Res>0.So, from (14) and (15), we have the following duality relations:
(16)Lφts=M0,mφts,1,Res>0,Mρ,mφts,τ=Lφτttm+τmρs,s,τ>0,Mρ,mφts,τ=1τLφttm/τm+τmρsτ,s,τ>0,(17)Mρ,mtmτm+τmρφts,τ=Lφτts,s,τ>0.Settingρ=0 in (14), we recover the natural transform defined as (see [49, 50])
(18)ℕφts,τ=∫0∞e−stφτtdt,s>0,τ>0.Thus, we have the followingMρ,m‐ℕ-transform duality
(19)ℕφts,τ=M0,mφts,τ,(20)Mρ,mφts,τ=ℕφttm/τm+τmρs,τ,s>0,τ>0,(21)Mρ,mtmτm+τmρφts,τ=ℕφts,τ;s>0,τ>0.The Sumudu transform is defined by [51–53]
(22)Sφtτ=∫0∞e−tφτtdt,τ>0.Thus,(23)Sφtτ=M0,mφt0,τ,τ>0,(24)Mρ,mφts,τ=1sSφttm/τm+τmρτs,s,τ>0.Based on these dualities of theMρ,m-transform (14) and these well-known integral transforms, it seems to be interesting to apply the Mρ,m-transform (14) in solving a variety of boundary and initial-boundary value problems. In this context, we recall the following results [48]:
(i)
Letφnt be the nth-order t-derivative of the function φt and φt≤Ket/γ with K>0, γ>0. Then,(25)Mρ,mtmτm+τmρφnts,τ=snτnℕφts,τ−∑k=0n−1skτk+1φn−k−10,where ℕφts,τ is defined by (18). Using the duality (21) in (25), we find
(26)ℕφnts,τ=snτnℕφts,τ−∑k=0n−1skτk+1φn−k−10,n=0,1,⋯.(ii)
Again, using the dualities stated before a convolution formula for theMρ,m-transform (14) can be obtained as follows. Here, the convolution for the Laplace transform will be considered; that is, for the functions φ and ψ, the convolution formula is given as(27)φ∗ψt=∫0tφxψt−xdx=∫0tφt−xψxdx.IfΦs,τ=Mρ,mφts,τ and Ψs,τ=Mρ,mψts,τ, thenτΦs,τΨs,τ=τ∫0∞e−st1φτt1/t1m+τmρdt1∫0∞e−st2ψτt2/t2m+τmρdt2=τ∫0∞∫0∞e−st1+t2φ~τt1ψ~τt2dt1dt2,where
(28)φ~t=φttm/τm+τmρ,ψ~t=ψttm/τm+τmρ.Settingt1+t2=t in the last equality, one gets
(29)τΦs,τΨs,τ=τ∫0∞∫t2∞e−stφ~τt−t2ψ~τt2dtdt2=τ∫0∞e−stdt∫0tφ~τt−t2ψ~τt2dt2.Here, changing of the integral order is used. Thus, using the duality of the Mρ,m and ℕ transforms (see (20)), we find
(30)τΦs,τΨs,τ=ℕφ~∗ψ~s,τ.Remark 2.
If we putρ=0 in (30), the case being interesting later in our work, then we get
(31)τℕφts,τℕψts,τ=ℕφ∗ψts,τ.(iii)
Once again, ...again, using the dualities stated before an inversion formula of theMρ,m-transform (14) is given (see Theorem 4.1 of [48]) as(32)φt=tmτm+τmρL−1Mρ,mφts,ttτ=12πitmτm+τmρ∫c−i∞c+i∞est/τMρ,mφts,τds,c,τ>0,as long as the integral converges absolutely. In case, when ρ=0 one obtains the following inversion formula of the natural transform (see Theorem 5.3 of [49])
(33)φt=12πi∫c−i∞c+i∞est/τℕφts,τds,c,τ>0.The residue theorem (see, e.g., [54]) is usually used to calculate the contour integrals in (32) and (33).


## 2. One-Layer Nonhomogeneous Diffusion System

Now, we investigate the solvability for the following one-layer nonhomogeneous initial boundary value problem:(34)∂φ∂t=d∂2φ∂x2+λt,τrx,t,x∈α,β,t,τ>0,(35)φx,0=ηx,x∈α,β,(36)ıφα,t+ι∂φ∂xα,t=λt,τζt,t≥0,τ>0,(37)ℓφβ,t+l∂φ∂xβ,t=λt,τξt,t≥0,τ>0,where d,ı,ι,ℓ, and l are constants such that ı+ι>0,ℓ+l>0, and λ,r,η,ζ, and ξ are given functions with λ as in (3).Applying theMρ,m-transform defined by (14), to (34), yields
(38)Mρ,m1λt,τφtx,t=dMρ,m1λt,τφxxx,t+Mρ,mrx,t.Using the duality of theMρ,m-transform and the natural transform given by (21) and (25), Equation (38) can be reduced to
(39)sτdℕφx,tx;s,τ−ℕφxxx,tx;s,τ−1τdφx,0=1dℕλt,τrx,tx;s,τ,where ℕφts,τ is defined by (18). Setting
(40)φ^x;s,τ=ℕφx,tx;s,τ,then, (39) can be expressed as
(41)φ^xxx;s,τ−sτdφ^x;s,τ=Fx;s,τ,where
(42)Fx;s,τ=−1dℕλt,τrx,tx;s,τ−1τdηx.Applying the variation of the parameter method to the nonhomogeneous equation (41) gives the general solution as
(43)φ^x;s,τ=Acoshsτdx+Bsinhsτdx+τds∫αxFy;s,τsinhsτdx−ydy,where A and B are arbitrary invariants which can depend on s and τ. Differentiating (43) with respect to x, gives
(44)φ^xx;s,τ=Asτdsinhsτdx+Bsτdcoshsτdx+∫αxFy;s,τcoshsτdx−ydy.Transforming the boundary conditions (36) and (37), implies
(45)ıφ^α;s,τ+ιφ^xα;s,τ=Mρ,mζt=ℕλt,τζt,ℓφ^β,,s,t+lφ^xβ;s,t=Mρ,mξt=ℕλt,τξt.For simplicity, we set the following vector notations:(46)a=ı,ι,b=ℓ,l,(47)Ly;s,τ=coshsτdy,sτdsinhsτdy,(48)Cy;s,τ=sinhsτdy,sτdcoshsτdy.Obviously, we have(49)∂L∂yy;s,τ=sτdCy;s,τ,∂C∂yy;s,τ=sτdLy;s,τ.Substituting (43) and (44) into (45) and using the vector notation, we give the algebraic linear system
(50)a,Lα;s,τa,Cα;s,τb,Lβ;s,τb,Cβ;s,τAB=Gs,τH+x;s,τ,where · is the usual dot product in ℝ2, and
(51)Gs,τ=ℕλt,τζt,(52)H+x;s,τ=Hs,τ−τds∫αβFy;s,τb,Cβ−y;s,τdy,with Hs,τ=ℕλt,τξt and ζt and ξt are the boundary data given in (36) and (37), respectively. The solution A,B of system (50) is
(53)AΔs=Gs,τb,Cβ;s,τ−H+x;s,τa,Cα;s,τ,BΔs=−Gs,τb,Lβ;s,τ−H+x;s,τa,Lα;s,τ,where
(54)Δs=b,Cβ;s,τa,Lα;s,τ−a,Cα;s,τb,Lβ;s,τ,is the determinant of the coefficient matrix of system (50). Substituting the constants A and B into (43) gives
(55)φ^x;s,τ=coshs/τdxΔsGs,τb,Cβ;s,τ−H+x;s,τa,Cα;s,τ+sinhs/τdxΔs−Gs,τb,Lβ;s,τ+H+x;s,τa,Lα;s,τ+τds∫αxsinhsτdx−yFy;s,τdy,which can be rewritten as
(56)φ^x;s,τ=Gs,τψx,β;s,τ,bΔs−Hs,τψx,α;s,τ,aΔs+θ^x;s,τ,where
(57)ψx,y;s,τ,L=L,Cy;s,τcoshsτdx−L,Ly;s,τsinhsτdx,(58)θ^x;s,τ=ψx,α;s,τ,aΔsτds∫αβb,Cβ−y;s,τFy;s,τdy+τds∫αxsinhsτdx−yFy;s,τdy.For further computation, we rewriteθ^x;s,τ as
(59)θ^x;s,τ=τds∫αxψx,α;s,τ,ab,Cβ−y;s,τ+Δssinhs/τdx−yΔsFy;s,τdy+τds∫xβψx,α;s,τ,ab,Cβ−y;s,τΔsFy;s,τdy.Lemma 3.
Letτ,s,x,y∈ℝ and L∈ℝ2. Then,
(60)ψx,y;s,τ,L=L,Cy−x;s,τ,(61)Δssinhsτdx−y=−ψx,α;s,τ,ab,Cβ−y;s,τ+a,Cα−y;s,τψx,β;s,τ,b.Consequently, for each zero s of the function Δssinhs/τdx−y, one has
(62)ψx,α;s,τ,ab,Cβ−y;s,τ=a,Cα−y;s,τψx,β;s,τ,b,where ψ is given by (57).Proof.
The first two conclusions of the lemma follow directly from the uniqueness theorem of the initial value problem for the second-order ordinary differential equations having constant coefficients.
For fixedτ,s,y∈ℝ and L∈ℝ2, in view of (46) and (57), the functions L,Cy−x;s,τ and ψx,y;s,τ,L are solutions to the following initial value problem:
(63)d2zxdx2−sτdzx=0,z0=L,Cy;s,τ,z′0=−sτdL,Ly;s,τ.
Thus, with the uniqueness of the solution to problem (63), we conclude (60).
It is easy to see that as functions inx, both sides of (61) are linear combinations of the functions coshs/τdx,andsinhs/τdx which are linearly independent solutions to the differential equation in (63). Thus, both sides of (61) solve this differential equation.
Moreover, in view of (57), both sides of (61) satisfy
(64)z0=−Δssinhsτdy,z′0=sτdcoshsτdy.
Hence, by the uniqueness theorem, (61) holds true.
For eachs∗ being a zero of the function Δssinhs/τdx−y, taking the limit in both sides of (61) as s⟶s∗ gives (62).
Applying Lemma3, (56) and (59) respectively can be reduced to
(65)φ^x;s,τ=Gs,τb,Cβ−x;s,τΔs−Hs,τa,Cα−x;s,τΔs+θ^x;s,τ,(66)θ^x;s,τ=τds∫αxa,Cα−y;s,τb,Cβ−x;s,τΔsFy;s,τdy+τds∫xβa,Cα−x;s,τb,Cβ−y;s,τΔsFy;s,τdy.
Next, in order to obtain the solution to the initial value problem (34)–(37), we apply the inversion formula (33) to (65) and (66). In doing so, we suppose that there are nonzero simple roots skk=1∞ of Δs. That is,
(67)Δsk=0,Δ′sk≠0,k=1,2,3,⋯.Lemma 4.
Suppose that (67) holds true. For each x,y∈ℝ and t,τ>0, we get
(68)ℕ−1a,Cx;s,τb,Cy;s,τs/τdΔs=Θx,y,t,τ,where
(69)Θx,y,t,τ=Θ0x,y,τ+∑k=1∞eskt/τa,Cx;sk,τb,Cy;sk,τsk/τdΔ′sk,with
(70)Θ0x,y=0ifa,Cx;s,τb,Cy;s,τs/τdΔs=O1,lims⟶0sτda,Cx;s,τb,Cy;s,τΔsifa,Cx;s,τb,Cy;s,τs/τdΔs=O1s.Proof.
Let(71)Rs,τ=a,Cx;s,τb,Cy;s,τs/τdΔs.
Applying the inversion formula (33), we find
(72)Θx,y,t=ℕ−1Rs,τ=12πi∫c−i∞c+i∞est/τRs,τds,c,τ>0.
The last integral can be usually calculated by the residue theorem [54]. Hence,
(73)Θx,y,t=∑polesskofRs,τResest/τRs,τ;sk.
Recalling (67), each skk=1,2,⋯ is a simple pole of est/τRs,τ. Therefore,
(74)Resest/τRs,τ;sk=eskt/τa,Cx;sk,τb,Cy;sk,τsk/τdΔ′sk,k=1,2,⋯.
Ats=0, we have
(75)Resest/τRs,τ;0=lims⟶0sest/τRs,τ=lims⟶0sΔsest/τa,Cx;s,τb,Cy;s,τs/τd=lims⟶0sτda,Cx;s,τb,Cy;s,τΔs.
We see that either(76)a,Cx;s,τb,Cy;s,τs/τdΔs=O1ora,Cx;s,τb,Cy;s,τs/τdΔs=O1s,as s tends to 0. Then, s=0 is either a removable singular point or a simple pole of est/τRs,τ.
Hence, substituting (74) and (75) in (73) gives the main conclusion of the Lemma, i.e., (69) and (70).☐In view of (68) and (42), we have
(77)ℕΘα−y,β−x,t=a,Cx−y,βb,Cβ−x,ss/τdΔs,where Θ is defined by (69) and (70), and
(78)Fy;s,τ=1dℕλy,ttm/τm+τmρy;s,τ−1τdηy.Hence, (66) can be rewritten as
(79)θ^x;s,τ=−1d∫αxℕΘα−y,β−x,tℕry,ttm/τm+τmρdy−1τd∫αxℕΘα−y,β−x,tηydy−1d∫xβℕΘα−x,β−y,tℕry,ttm/τm+τmρdy−1τd∫xβℕΘα−x,β−y,tηydy.By the convolution formula (31), the inverse natural transform of (79) is
(80)θx,t,τ=−1τd∫αx∫0tΘα−y,β−x,t−ζry,ζζm/τm+τmρdζdy−1τd∫αxΘα−y,β−x,tηydy−1τd∫xβ∫0tΘα−x,β−y,t−ζry,ζζm/τm+τmρdζdy−1τd∫xβΘα−x,β−y,tηydy.That is,(81)θx,t,τ=−1τd∫αxΘ0α−y,β−xηydy−1τd∫xβΘ0α−x,β−yηydy−1τd∫αx∫0tΘ0α−y,β−xry,ζζm/τm+τmρdζdy−1τd∫xβ∫0tΘ0α−x,β−yry,ζζm/τm+τmρdζdy−1τd∑k=1∞b,Cβ−x;sk,τsk/τdΔ′sk∫αx∫0teskt−ζ/τry,ζa,Cα−y;sk,τζm/τm+τmρdζdy−1τd∑k=1∞a,Cα−x;sk,τsk/τdΔ′sk∫xβ∫0teskt−ζ/τry,ζb,Cβ−y;sk,τζm/τm+τmρdζdy−1τd∑k=1∞eskt/τb,Cβ−x;sk,τsk/τdΔ′sk∫αxa,Cα−y;sk,τηydy−1τd∑k=1∞eskt/τa,Cα−x;sk,τsk/τdΔ′sk∫xβb,Cβ−y;sk,τηydy.From Lemma3, one has
(82)a,Cα−x;s,τb,Cβ−y;s,τ=a,Cα−y;s,τb,Cβ−x;s,τ,at s=0 and s=skk=1,2,⋯ the zeros of Δssinhs/τdx−y. That results in
(83)Θ0α−y,β−x=Θ0α−x,β−y,a,Cα−x;sk,τ∫xβb,Cβ−y;sk,τϑydy=b,Cβ−x;sk,τ∫xβa,Cα−y;sk,τϑydy.The first conclusion is obvious whenΘ0=0 in (70). Thus, (81) can be simplified as
(84)θx,t,τ=−1τd∫αβΘ0α−x,β−yηydy−1τd∫αβ∫0tΘ0α−x,β−yry,ζζm/τm+τmρdζdy−1τd∑k=1∞b,Cβ−x;sk,τsk/τdΔ′sk∫αβ∫0teskt−ζ/τry,ζa,Cα−y;sk,τζm/τm+τmρdζdy−1τd∑k=1∞eskt/τb,Cβ−x;sk,τsk/τdΔ′sk∫αβa,Cα−y;sk,τηydy.Next, we return to (65). Using (26) (for n=1) and (51), (65) can be rewritten as
(85)φ^x;s,τ=sℕλt,τζtb,Cβ−x;s,τsΔs−sℕλt,τξta,Cα−x;s,τsΔs+θ^x;s,τ=τℕddtλt,τζt+λ0,τζ0b,Cβ−x;s,τsΔs−τℕddtλt,τξt+λ0,τξ0a,Cα−x;s,τsΔs+θ^x;s,τ,where θ^x;s,τ is given in (66). Now, we can obtain the solution φx,t of Problem (34)–(37) by operating the inversion formula (33) in (85). In doing so, we need the following lemma.Lemma 5.
Assume that (67) holds true. Then, for each y∈ℝ,t,τ>0 and L∈ℝ2, we get
(86)ℕ−1L,Cy;s,τsΔs=Φy,t;L,where
(87)Φy,t;L=Φ0y;L+∑k=1∞eskt/τL,Cy;sk,τskΔ′sk,(88)Φ0y;L=lims⟶0L,Cy;s,τΔs,ifL,Cy;s,τsΔs=O1slims⟶0∂∂ssL,Cy;s,τΔs,ifL,Cy;s,τsΔs=O1s2.Proof.
The proof is similar to Lemma4.
From Lemma5, we see thatℕ−1b,Cβ−x;s,τ/sΔs=Φβ−x,t;b,ℕ−1a,Cα−x;s,τ/sΔs=Φα−x,t;a.
Hence, in view of the convolution formula (31) and the inversion of natural transform (33), inverting (85) yields
(89)φx,t,τ=∫0tΦβ−x,t−ς;bζ~′ς+ζ~0δ0ςdς−∫0tΦα−x,t−ς;aξ~′ς+ξ~0δ0ςdς+θx,t,τ,where ζ~=λt,τζt,ξ~=λt,τξt,δ0 is the well-known Dirac delta function, and θx,t,τ is given by (84). Then, using the basic property of the Dirac delta function, that is, δ0ςΦς=Φ0, results in
(90)φx,t,τ=∫0tΦβ−x,t−ς;bζ′~ςdς+ζ~0Φβ−x,t;b−∫0tΦα−x,t−ς;aξ′~ςdς−ξ~0Φα−x,t;a+θx,t,τ.
Integrating by parts gives(91)φx,t,τ=λt,τζtΦβ−x,0;b−∫0tλς,τD2Φβ−x,t−ς;bζςdς−λt,τΦα−x,0;aξt+∫0tλς,τD2Φα−x,t−ς;aξςdς+θx,t,τ,where D2=∂/∂tΦx,t;b. Substituting from (87) gives
(92)φx,t,τ=λt,τζtΦ0β−x;b+∑k=1∞b,Cβ−x;sk,τskΔ′sk−∫0tλς,τζς∑k=1∞eskt−ς/τb,Cβ−x;sk,τΔ′skdς−λt,τξtΦ0α−x;a+∑k=1∞a,Cα−x;sk,τskΔ′sk+∫0tλς,τξς∑k=1∞eskt−ς/τa,Cα−x;sk,τΔ′skdς+θx,t,τ,with θx,t,τ is given by (84). This result can be rewritten as
(93)φx,t,τ=λt,τζtΦ0β−x;b+∑k=1∞ΓkζtskΔ′skb,Cβ−x;sk,τ−λt,τξtΦ0α−x;a−∑k=1∞ΓkξtskΔ′ska,Cα−x;sk,τ+θx,t,τ,where Γk is the operator defined as
(94)Γkϕt=λt,τϕt−sk∫0tλς,τeskt−ς/τϕςdς.The integral in (94) is the Laplacian convolution formula for λt,τϕt with eskt/τ. As a result, (93), together with (84) and (94), expresses the solution of Problem (34)–(37).Remark 6.
Whenρ=0 and rj=νj=0, for all j=1,⋯,n, Problem (34)–(37) and its solution
(95)φx,t,τ=1=ζtΦ0β−x;b+∑k=1∞Γ~kζtskΔ′skb,Cβ−x,sk−ξtΦ0α−x;a−∑k=1∞Γ~kξtskΔ′ska,Cα−x,sk+θx,t,with Φ0α−x;a,Cy,sk=Cy;sk,τ=1 defined as (88), (46), respectively,
(96)Γ~kφt=φt−sk∫0teskt−ςφςdς,θx,t=−1d∫αβΘ0α−x,β−yηydy−∑k=1∞esktb,Cβ−x;sk,1dskΔ′sk∫αβa,Cα−y;sk,1ηydy,are reduced to that in Section 3 of [38].


### 2.1. Illustrative Examples

Here, we discuss two illustrative test cases to show the accuracy and effectiveness of our technique.Example 1.
Heat equation with zero temperatures at finite ends.
The following initial boundary value problem with homogeneous Dirichlet boundary conditions(97)∂φ∂t=d∂2φ∂x2,x∈0,L,t>0,φx,0=fx,x∈0,L,φ0,t=0,φL,t=0,t≥0,is a special case of the one-layer diffusion system (34)–(37) when α=0,β=L,ı=1,ι=0,ℓ=1,l=0 and rx,t=ζt=ξt=0,ηx=fx. For simplicity, we will take τ=1. Thus, from (93), we have φx,t=θx,t. From (46), we have a=b=1,0,
(98)Lx;s=coshsdx,sdsinhsdx,Cx;s=sinhsdx,sdcoshsdx.
Hence, (54) yields
(99)Δs=sinhsdL,Δ′s=L2ddscoshsdL.
So, we have(100)Δsk=0,Δ′sk≠0,k=1,2,3,⋯,at
(101)sk=−k2π2dL2,k=1,2,3,⋯.
Moreover,(102)a,Cx;sb,Cy;ss/dΔs=sinhs/dxsinhs/dys/dsinhs/dL=O1,as s⟶0. Therefore, (70) gives Θ0x,y=0. Further,
(103)a,C−y;sk=−isinkπLy,b,CL−x;sk=−1k+1isinkπLx,k=1,2,3,⋯,skdΔ′sk=−1kL2d,k=1,2,3,⋯.
Finally, from (80), we get
φx,t=θx,t=−1/d∑k=1∞esktb,CL−x;sk/sk/dΔ′sk∫0La,C−y;skfydy=∑k=0∞2/L∫0Lfysinkπy/Ldye−k2π2d/L2tsinkπ/Lx,which recovers the solution to problem (97) obtained via the separation of variables method in [55].Example 2.
Heat flow with sources and homogeneous boundary conditions.
The following initial boundary value problem(104)∂φ∂t=d∂2φ∂x2+tm+1−ρe−tsin3x,x∈0,π,t>0,φx,0=fx,x∈0,π,φ0,t=0.φπ,t=0,t≥0,is a special case of the one-layer diffusion system (34)–(37) when α=0,β=π,d=τ=1,ı=1,ι=0,ℓ=1,l=0 and rx,t=e−tsin3x,ηx=fx,ζt=ξt=0. Thus, from (93), we have φx,t=θx,t. From (46), we have a=b=1,0,
(105)Lx;s=coshsx,ssinhsx,Cx;s=sinhsx,scoshsx.
Hence, (54) yields
(106)Δs=sinhsπ,Δ′s=π21scoshsπ.
So,∆s has simple zeros at sk=−k2,k=1,2,3,⋯. Further,
(107)a,C−y;sk=−isinky,b,Cπ−x;sk=−1k+1isinkx,k=1,2,3,⋯,skΔ′sk=−1kπ2,k=1,2,3,⋯.
On the other hand,(108)a,Cx;sb,Cy;ssΔs=sinhsxsinhsyssinhsπ=O1as s⟶0. Therefore, (70) gives Θ0x,y=0.
Finally, from (80), we get
(109)φx,t=θx,t=∑k=1∞2π∫0π∫0tek2−1zsin3ysinkyzm+1ρdzdye−k2tsinkx+∑k=0∞2π∫0πfysinkydye−k2tsinkx=∫0te8zzm+1ρe−9tsin3x+∑k=0∞2π∫0πfysinkydye−k2tsinkx=∑k=o∞aktsinkx,where
(110)akt=ak0e−k2t,k≠3,a30+ate−9t,k=3,with
(111)ak0=2π∫0tfysinkydy,at=2π∫0te8zzm+1ρdz.Whenρ=0, formula (109) recovers the solution to problem (104) (when ρ=0) obtained via the eigenfunction expansion method in [55].


## 2.1. Illustrative Examples

Here, we discuss two illustrative test cases to show the accuracy and effectiveness of our technique.Example 1.
Heat equation with zero temperatures at finite ends.
The following initial boundary value problem with homogeneous Dirichlet boundary conditions(97)∂φ∂t=d∂2φ∂x2,x∈0,L,t>0,φx,0=fx,x∈0,L,φ0,t=0,φL,t=0,t≥0,is a special case of the one-layer diffusion system (34)–(37) when α=0,β=L,ı=1,ι=0,ℓ=1,l=0 and rx,t=ζt=ξt=0,ηx=fx. For simplicity, we will take τ=1. Thus, from (93), we have φx,t=θx,t. From (46), we have a=b=1,0,
(98)Lx;s=coshsdx,sdsinhsdx,Cx;s=sinhsdx,sdcoshsdx.
Hence, (54) yields
(99)Δs=sinhsdL,Δ′s=L2ddscoshsdL.
So, we have(100)Δsk=0,Δ′sk≠0,k=1,2,3,⋯,at
(101)sk=−k2π2dL2,k=1,2,3,⋯.
Moreover,(102)a,Cx;sb,Cy;ss/dΔs=sinhs/dxsinhs/dys/dsinhs/dL=O1,as s⟶0. Therefore, (70) gives Θ0x,y=0. Further,
(103)a,C−y;sk=−isinkπLy,b,CL−x;sk=−1k+1isinkπLx,k=1,2,3,⋯,skdΔ′sk=−1kL2d,k=1,2,3,⋯.
Finally, from (80), we get
φx,t=θx,t=−1/d∑k=1∞esktb,CL−x;sk/sk/dΔ′sk∫0La,C−y;skfydy=∑k=0∞2/L∫0Lfysinkπy/Ldye−k2π2d/L2tsinkπ/Lx,which recovers the solution to problem (97) obtained via the separation of variables method in [55].Example 2.
Heat flow with sources and homogeneous boundary conditions.
The following initial boundary value problem(104)∂φ∂t=d∂2φ∂x2+tm+1−ρe−tsin3x,x∈0,π,t>0,φx,0=fx,x∈0,π,φ0,t=0.φπ,t=0,t≥0,is a special case of the one-layer diffusion system (34)–(37) when α=0,β=π,d=τ=1,ı=1,ι=0,ℓ=1,l=0 and rx,t=e−tsin3x,ηx=fx,ζt=ξt=0. Thus, from (93), we have φx,t=θx,t. From (46), we have a=b=1,0,
(105)Lx;s=coshsx,ssinhsx,Cx;s=sinhsx,scoshsx.
Hence, (54) yields
(106)Δs=sinhsπ,Δ′s=π21scoshsπ.
So,∆s has simple zeros at sk=−k2,k=1,2,3,⋯. Further,
(107)a,C−y;sk=−isinky,b,Cπ−x;sk=−1k+1isinkx,k=1,2,3,⋯,skΔ′sk=−1kπ2,k=1,2,3,⋯.
On the other hand,(108)a,Cx;sb,Cy;ssΔs=sinhsxsinhsyssinhsπ=O1as s⟶0. Therefore, (70) gives Θ0x,y=0.
Finally, from (80), we get
(109)φx,t=θx,t=∑k=1∞2π∫0π∫0tek2−1zsin3ysinkyzm+1ρdzdye−k2tsinkx+∑k=0∞2π∫0πfysinkydye−k2tsinkx=∫0te8zzm+1ρe−9tsin3x+∑k=0∞2π∫0πfysinkydye−k2tsinkx=∑k=o∞aktsinkx,where
(110)akt=ak0e−k2t,k≠3,a30+ate−9t,k=3,with
(111)ak0=2π∫0tfysinkydy,at=2π∫0te8zzm+1ρdz.Whenρ=0, formula (109) recovers the solution to problem (104) (when ρ=0) obtained via the eigenfunction expansion method in [55].


## 3. Multilayer Nonhomogeneous Diffusion System

Here, we are seeking the solution of our main problem defined in (2)–(8), which was converted into a sequence of initial boundary value problems (9)–(11). For the convenience of the reader and in order to draw the full picture in an easy way, we start with solving the bilayer diffusion problem in the following subsection; then, we move to the general case in Section 3.2.


### 3.1. Solution of a Two-Layer Problem

For the two-layer problem, we have(112)∂φ1∂t=d1∂2φ1∂x2+λt,τr1x,t,x∈x0,x1,t,τ>0,(113)φ1x,0=η1x,x∈x0,x1,(114)ıφ1x0,t+ι∂φ1∂xx0,t=λt,τζ1t,t≥0,τ>0,(115)v1φ1x1,t+μ1∂φ1∂xx1,t=λt,τξ1t,t≥0,τ>0,(116)∂φ2∂t=d2∂2φ2∂x2+λt,τr2x,t,x∈x1,x2,t,τ>0,(117)φ2x,0=η2x,x∈x1,x2,(118)v2φ2x1,t+μ2∂φ2∂xx1,t=λt,τζ2t,t≥0,τ>0,(119)ℓφ2x2,t+l∂φ2∂xx2,t=λt,τξ2t,t≥0,τ>0.Similar to what we denote in Section2, we define the following vector notation a1=ı,ι,b1=ν1,µ1,a2=ν2,µ2,b2=ℓ,l, and
(120)L1y;s,τ=coshsτd1y,sτd1sinhsτd1y,L2y;s,τ=coshsτd2y,sτd2sinhsτd2y,C1y;s,τ=sinhsτd1y,sτd1coshsτd1y,C2y;s,τ=sinhsτd2y,sτd2coshsτd2y.Also, analogues to (54), define
(121)Δ1s=b1,C1x1,s,τa1,L1x0,s,τ−a1,C1x0,s,τb1,L1x1,s,τ,Δ2s=b2,C2x2,s,τa2,L2x1,s,τ−a2,C2x1,s,τb2,L2x2,s,τ.Further, similar to (67), suppose that there are nonzero simple roots sk1k=1∞ and sk2k=1∞ of the functions Δ1sandΔ2s, respectively. That is,
(122)Δ1sk1=0,Δ1′sk1≠0,Δ2sk2=0,Δ2′sk2≠0k=1,2,⋯.Therefore, according to (84), we obtain
(123)θ1x,t,τ=−1τd1∫x0x1Θ01x0−x,x1−yη1ydy−1τd1∫x0x1∫0tΘ01x0−x,x1−yr1y,ςςm/τm+τmρdςdy−1τd1∑k=1∞b1,C1x1−x,sk1,τsk1/τd1Δ1′sk1∫x0x1∫0tesk1t−ς/τr1y,ςa1,C1x0−y,sk1,τςm/τm+τmρdςdy−1τd1∑k=1∞esk1t/τb1,C1x1−x,sk1,τsk1/τd1Δ1′sk1∫x0x1a1,C1x0−y,sk1,τη1ydy,(124)θ2x,t,τ=−1τd2∫x1x2Θ02x1−x,x2−yη2ydy−1τd2∫x1x2∫0tΘ02x1−x,x2−yr2y,ςςm/τm+τmρdςdy−1τd2∑k=1∞b2,C2x2−x,sk2,τsk2/τd2Δ2′sk2∫x1x2∫0tesk2t−ς/τr2y,ςa2,C2x1−y,sk2,τςm/τm+τmρdςdy−1τd2∑k=1∞esk2t/τb2,C2x2−x,sk2,τsk2/τd2Δ2′sk2∫x1x2a2,C2x1−y,sk2,τη2ydy,whereΘ01andΘ02 can be defined as in Lemma 4.Also, similar to (93), with the respective forms Φ01andΦ02 from Lemma 5 and the matching condition ξ1t=ζ2t, we get
(125)φ1x,t,τ=λt,τζ1tΦ01x1−x;b1+∑k=1∞Γk1ζ1tsk1Δ1′sk1b1,C1x1−x,sk1,τ−λt,τξ1tΦ01x0−x;a1−∑k=1∞Γk1ξ1tsk1Δ1′sk1a1,C1x0−x,sk1,τ+θ1x,t,τ,(126)φ2x,t,τ=λt,τξ1tΦ02x2−x;b2+∑k=1∞Γk2ξ1tsk2Δ2′sk2b2,C2x2−x,sk2,τ−λt,τξ2tΦ02x1−x;a2−∑k=1∞Γk2ξ2tsk2Δ2′sk2a2,C2x1−x,sk2,τ+θ2x,t,τ,where the operators Γk1andΓk2 are obtained from (94). The matching condition φ1x1,t=Λ1φ2x1,t yields
(127)Λ1λt,τΦ02x2−x1;b2+λt,τΦ01x0−x1;a1ξ1t+∑k=1∞a1,C1x0−x1,sk1,τsk1Δ1′sk1Γk1ξ1t−Λ1b2,C2x2−x1,sk2,τsk2Δ2′sk2Γk2ξ1t=λt,τζ1tΦ010;b1+θ1x1,t,τ+∑k=1∞b1,C10,sk1,τsk1Δ1′sk1Γk1ζ1t+Λ1λt,τξ2tΦ020;a2+∑k=1∞Λ1a2,C20,sk2,τsk2Δ2′sk2Γk2ξ2t−Λ1θ2x1,t,τ.For the unknown functionξ1, we can rewrite the linear integral equation (127) as
(128)λt,τξ1t+∑k=1∞akΓk1ξ1t+bkΓk2ξ1t=ct,where
(129)ak=a1,C1x0−x1,sk1,τsk1Δ1′sk1Λ1Φ02x2−x1;b2+Φ01x0−x1;a1,bk=−Λ1b2,C2x2−x1,sk2,τsk2Δ2′sk2Λ1Φ02x2−x1;b2+Φ01x0−x1;a1,ct=1Λ1Φ02x2−x1;b2+Φ01x0−x1;a1λt,τζ1tΦ010;b1+θ1x1,t,τ+∑k=1∞b1,C10,sk1,τsk1Δ1′sk1Γk1ζ1t+Λ1λt,τξ2tΦ020;a2+∑k=1∞Λ1a2,C20,sk2,τsk2Δ2′sk2Γk2ξ2t−Λ1θ2x1,t,τ.Inspired by the convolution formula (31), the natural transform of (128) is
(130)ℕξ1tλt,τs,τ+∑k=1∞akℕξ1tλt,τs,τ−τsk1s−sk1ℕξ1tλt,τ+bkℕξ1tλt,τs,τ−τsk2s−sk2ℕξ1tλt,τ=ℕcts,τ,which can be rewritten as
(131)1−ℕψts,τℕξ1tλt,τs,τ=ℕcts,τ.That is,(132)ℕξ1tλt,τs,τ=11−ℕψts,τℕcts,τ=1+ℕψts,τ+ℕψts,τ2+⋯ℕcts,τ=ℕcts,τ+ℕcts,τ∑m=1∞ℕψts,τm,where
(133)ℕψts,τ=−∑k=1∞ak1−τsk1s−sk1+bk1−τsk2s−sk2,for which the inverse natural transform is
(134)ψt=−τ∑k=1∞ak+bkδ0t+τ∑k=1∞aksk1esk1/τt+bksk2esk2/τt,where δ0 is the Dirac delta function. Hence, we have
(135)ξ1t=1λt,τct+∑m=1∞ψm∗ct,where for m≥2,ψm is the m-times self-convolution of ψ. Thus, one can conclude the solution to the bilayer diffusion problem (112)–(119) by the formulas (125) and (126), together with (123) and (124), with ξ1 given in (135). Now, it is time to attack an illustrative example in the following subsection.


#### 3.1.1. Illustrative Example

Here, we discuss the solvability of the following two-layer diffusion system.Example 3.
Temperature distribution in the two-layer slab with mixed boundary condition.
Consider the following initial boundary value system. The diffusion equations are(136)∂φk∂t=dk∂2φk∂x2,xk−1<x<xk,x0=0,t>0,k=1,2,the initial conditions are
(137)φkx,0=ηkx,xk−1≤x≤xk,k=1,2,the outer boundary conditions are
(138)φ10,t=0,∂φ2x2,t∂x=0,t≥0,and the interface conditions at x=x1 are
(139)φ1x1,t=φ2x1,t,d1∂φ1x1,t∂x=d2∂φ2x1,t∂x,t≥0.
Comparing this problem with the general one that is defined by (112)–(119) reveals
(140)a1=1,0,b1=0,d1,a2=0,d2,b2=0,1,r1x,t=r2x,t=ζ1t=ξ2t=0,ξ1t=ζ2t.Here, we will consider ρ=0andτ=1.
Then, the solution of the problem (136)–(139) can be obtained from (125) and (126) as
(141)φ1x,t=−ξ1tΦ01x0−x;a1−∑k=1∞Γk1ξ1tsk1Δ1′sk1a1,C1x0−x,sk1,τ+θ1x,t,(142)φ2x,t=ξ1tΦ02x2−x;b2+∑k=1∞Γk2ξ1tsk2Δ2′sk2b2,C2x2−x,sk2,τ+θ2x,t,with
(143)θ1x,t=−1d1∫x0x1Θ01x0−x,x1−yη1ydy−1d1∑k=1∞esk1tb1,C1x1−x,sk1sk1/d1Δ1′sk1∫x0x1a1,C1x0−y,sk1η1ydy,θ2x,t=−1d2∫x1x2Θ02x1−x,x2−yη2ydy−1d2∑k=1∞esk2tb2,C2x2−x,sk2sk2/d2Δ2′sk2∫x1x2a2,C2x1−y,sk2η2ydy.The operatorsΓk1andΓk2 are obtained from (94).Next, we are going to simplify these formulas. Direct computations give(144)Δ1s=d1sd1coshsd1x1,Δ2s=−d2sd2sinhsd2x2−x1,Δ1′s=d1x12d1sinhsd1x1+12d1scoshsd1x1,Δ2′s=−d21d2sinhsd2x2−x1+x2−x12d2sd2coshsd2x2−x1.It is clear thatΔ1sk1=0andΔ2sk2=0 when
(145)sk1=2k−12π2d14x12,sk2=−k2π2d2x2−x12,k=1,2,3,⋯,respectively. Moreover,
(146)a1,C1x;sb1,C1y;ss/d1Δ1s=sinhs/d1xcoshs/d1ys/d1coshs/d1x1=O1,a2,C2x;sb2,C2y;ss/d2Δ2s=−coshs/d2xcoshs/d2ys/d2sinhs/d2x2−x1=O1s,as s⟶0. Therefore, (70) gives Θ01x,y=0andΘ02x,y=−d2/x2−x1.Further,(147)a1,C1−y;sk1=−isin2k−1π2x1y,b1,C1x1−x;sk1=−1k+12k−1πd12x1isin2k−1π2x1x,sk1d1Δ1′sk1=−1k2k−1π4,k=1,2,3,⋯,a2,C2x1−y;sk2=kπd2ix2−x1coskπx2−x1x1−y,b2,C2x2−x;sk2=kπix2−x1coskπx2−x1x2−x,sk2d2Δ2′sk2=−1kk2π22x2−x1.Thus,(148)θ1x,t=∑k=1∞2x1∫0x1η1ysin2k−1πy2x1dye−2k−12π2/4x12d1tsin2k−1π2x1x,θ2x,t=−1x2−x1∫x1x2η2ydy+∑k=1∞−1k1x2−x1∫x1x2η2ycoskπx2−x1x1−ydye−k2π22π2/x2−x12d2tcoskπx2−x1x2−x1.Similarly, we have(149)a1,C1x0−x,ssΔ1s=sinhs/d1x0−xsd1s/d1coshs/d1x1=O1s,b2,C2x2−x,ssΔ2s=−s/d2coshs/d2x2−xsd2s/d1sinhs/d2x2−x1=O1s2.So, (88) gives
(150)Φ01x0−x;a1=−xd1,Φ02x2−x;b2=16d2x2−x13x2−x2−x2−x12.Hence, loading the quantities in (141) and (142) gives the solution to problem (136)–(139) as
(151)φ1x,t=xd1ξ1t+8x1π2d1∑k=1∞−1kΓk1ξ1t2k−12sin2k−1π2x1x+θ1x,t,φ2x,t=−3x2−x2−x2−x126d2x2−x1ξ1t+2x2−xd2π2∑k=1∞−1kΓk2ξ1tk2coskπx2−x1x2−x+θ2x,t.Applying the matching conditionφ1x1,t=φ2x1,t gives
(152)ξ1t+∑k=1∞akΓk1ξ1t+bkΓk2ξ1t=ct,where
(153)ak=−24d2x12k−12π23d2−d1x1+d1x2,bk=−6d1x2−x1k2π23d2−d1x1+d1x2,ct=3d1d23d2−d1x1+d1x2θ1x1,t−θ2x1,t.To solve this integral equation, we use (134), to obtain
(154)ψt=δ0t−6d1d23d2−d1x1+d1x2∑k=1∞1x1e−2k−12π2d1/4x12t+1x2−x1e−k2π2d2/x2−x12t,where δ0 is the Dirac delta function. Here, we used
(155)∑k=1∞12k−12=π28,∑k=1∞1k2=π26.Hence, from (135), we have
(156)ξ1t=ct+ψ∗ct+ψ∗ψ∗ct+⋯=ct+∑m=1∞ψm∗ct.


### 3.2. Solution of a Multilayer Problem

Here, we investigate the solvability of the main problem (2)–(8), through solving the initial boundary value problems (9)–(11). Similar to what we have denoted in Section 2, we consider the following notations:
(157)aj=ı,ι,j=1,vj,μj,2≤j≤n,bj=vj,μj,1≤j≤n−1,ℓ,l,j=n,and, for all 1≤j≤n,
(158)Ljy;s,τ=coshsτdjy,sτdjsinhsτdjy,Cjy;s,τ=sinhsτdjy,sτdjcoshsτdjy.Moreover, we define(159)Δjs=bj,Cjxj,s,τaj,Ljxj−1,s,τ−aj,Cjxj−1,s,τbj,Ljxj,s,τ,and let skjk=1∞ be the sequence of zeros of the function Δjsforall1≤j≤n, i.e.,
(160)Δjskj=0,Δj′skj≠0k=1,2,⋯.Analogue to the computations of (84) and (93), we have for the current case, for all j=1,⋯,n,
(161)θjx,t,τ=−1τdj∫xj−1xjΘ0jxj−1−x,xj−yηjydy−1τdj∫xj−1xj∫0tΘ0jxj−1−x,xj−yrjy,ζζm/τm+τmρdζdy−1τdj∑k=1∞bj,Cjxj−x,skj,τskj/τdjΔj′skj∫xj−1xj∫0teskjt−ζ/τrjy,ζaj,Cjxj−1−y,skj,τζm/τm+τmρdζdy−1τdj∑k=1∞eskjt/τbj,Cjxj−x,skj,τskj/τdjΔj′skj∫xj−1xjaj,Cjxj−1−y,skj,τηjydy,where Θ0j can be defined in a similar way as in Lemma 4, and
(162)φjx,t,τ=λt,τζjtΦjxj−x,0;bj−∫0tλς,τD2Φjxj−x,t−ς;bjζjςdς−λt,τξjtΦjxj−1−x,0;aj+∫0tλς,τD2Φjxj−1−x,t−ς;ajξjςdς+θjx,t,τ,with the respective forms Φj defined by (87) in Lemma 5. This last equation (162) can be rewritten as
(163)φjx,t,τ=Tjζ~jxj−x,t;bj−Tjξ~jxj−1−x,t;aj+θjx,t,τ,in which ζ~j=λt,τζjt,ξ~j=λt,τξjt and the linear operator Tj is defined by
(164)Tjφy,t;L=Φjy,0;Lφt−∫0tD2Φjy,t−ς;Lφςdς,for all j=1,⋯,n,L∈ℝ2.The matching conditionsφjxj,t,τ=Λjφj+1xj,t,τ,j=1,⋯,n−1, lead to
(165)Tjζ~j0,t;bj−Tjξ~jxj−1−xj,t;aj+θjxj,t,τ=ΛjTj+1ζ~j+1xj+1−xj,t;bj+1−Tj+1ξ~j+10,t;aj+1+θj+1xj,t,τ.Using the matching conditions (12), we have ζ~j+1t=ξ~jt for all 1≤j≤n−1. Thus, for j=1,
(166)−T1ξ~1x0−x1,t;a1−Λ1T2ξ~1x2−x1,t;b2+Λ1T2ξ~20,t;a2=Λ1θ2x1,t,τ−θ1x1,t,τ−T1ζ~10,t;b1.For2≤j≤n−2,
(167)Tjξ~j−10,t;bj−Tjξ~jxj‐1−xj,t;aj−ΛjTj+1ξ~jxj+1−xj,t;bj+1+ΛjTj+1ξ~j+10,t;aj+1=Λjθj+1xj,t,τ−θjxj,t,τ.Forj=n−1,
(168)Tn−1ξ~n−20,t;bn−1−Tn−1ξ~n−1xn−2−xn−1,t;an−1−Λn−1Tnξ~n−1xn−xn−1,t;bn=Λn−1θnxn−1,t,τ−θn−1xn−1,t,τ−Λn−1Tnξ~n0,t;an.System (166), (167), and (168), of n−1 integral equations of the unknowns ξ~j;1≤j≤n−1, can be adjusted as a matrix equation
(169)A0ht+A′∗ht=bt,with At is a tridiagonal matrix of order n−1 whose entries are as follows:
(170)−Φjxj−1−xj,t;aj−ΛjΦjxj+1−xj,t;bj+1,1⩽j⩽n−1maindiagonal,ΛjΦj+10,t;aj+1,1⩽j⩽n−2superdiagonalΦj0,t;bj,2⩽j⩽n−1subdiagonal,and the vectors ht and bt are defined as
(171)ht=ξ~1t⋮ξ~n−1t,bt=Λ1θ2x1,t,τ−θ1x1,t,τ−T1ζ~10,t;b1Λ2θ3x2,t,τ−θ2x2,t,τ⋮Λn−2θn−1xn−2,t,τ−θn−2xn−2,t,τΛn−1θnxn−1,t,τ−θn−1xn−1,t,τ−Λn−1Tnξ~n0,t;an.In fact, we can rewrite (169) as
(172)ht=Ct+B∗ht,with Ct=A0−1bt and Bt=−A0−1A′t. In view of the convolution formula (31), the natural transform of (172) reads
(173)ℕhts,τ=ℕCts,τ+τℕBts,τℕhts,τ,which is equivalent to
(174)ℕhts,τ=I−τℕBts,τ−1ℕCts,τ=I+τℕBts,τ+τℕBts,τ2+⋯ℕCts,τ,where I is the n−1×n−1 identity matrix. Once again, throughout the convolution sense (31), the natural transform inversion of (174) is
(175)ht=Ct+B∗Ct+B∗B∗Ct+⋯=Ct+∑m=1∞Bm∗Ct,where Bm is the m-times self-convolution of B. Finally, the solution of the nonhomogeneous multilayer diffusion systems (9)–(11) and hence that of the main problem (2)–(8) is concluded as
(176)φjx,t,τ=λt,τζjtΦ0jxj−x;bj−λt,τξjtΦ0jxj−1−x;aj+∑k=1∞λt,τζjt−∫0tλς,τskjeskjt−ς/τζjςdςskj∆j′skjbj,Cjxj−x,skj,τ−∑k=1∞λt,τξjt−∫0tλς,τskjeskjt−ς/τξjςdςskj∆j′skjaj,Cjxj−1−x,skj,τ+θjx,t,τ,with the respective forms Φ0j and θj defined as in (88) and (161), respectively, for all j=1,⋯,n.


## 3.1. Solution of a Two-Layer Problem

For the two-layer problem, we have(112)∂φ1∂t=d1∂2φ1∂x2+λt,τr1x,t,x∈x0,x1,t,τ>0,(113)φ1x,0=η1x,x∈x0,x1,(114)ıφ1x0,t+ι∂φ1∂xx0,t=λt,τζ1t,t≥0,τ>0,(115)v1φ1x1,t+μ1∂φ1∂xx1,t=λt,τξ1t,t≥0,τ>0,(116)∂φ2∂t=d2∂2φ2∂x2+λt,τr2x,t,x∈x1,x2,t,τ>0,(117)φ2x,0=η2x,x∈x1,x2,(118)v2φ2x1,t+μ2∂φ2∂xx1,t=λt,τζ2t,t≥0,τ>0,(119)ℓφ2x2,t+l∂φ2∂xx2,t=λt,τξ2t,t≥0,τ>0.Similar to what we denote in Section2, we define the following vector notation a1=ı,ι,b1=ν1,µ1,a2=ν2,µ2,b2=ℓ,l, and
(120)L1y;s,τ=coshsτd1y,sτd1sinhsτd1y,L2y;s,τ=coshsτd2y,sτd2sinhsτd2y,C1y;s,τ=sinhsτd1y,sτd1coshsτd1y,C2y;s,τ=sinhsτd2y,sτd2coshsτd2y.Also, analogues to (54), define
(121)Δ1s=b1,C1x1,s,τa1,L1x0,s,τ−a1,C1x0,s,τb1,L1x1,s,τ,Δ2s=b2,C2x2,s,τa2,L2x1,s,τ−a2,C2x1,s,τb2,L2x2,s,τ.Further, similar to (67), suppose that there are nonzero simple roots sk1k=1∞ and sk2k=1∞ of the functions Δ1sandΔ2s, respectively. That is,
(122)Δ1sk1=0,Δ1′sk1≠0,Δ2sk2=0,Δ2′sk2≠0k=1,2,⋯.Therefore, according to (84), we obtain
(123)θ1x,t,τ=−1τd1∫x0x1Θ01x0−x,x1−yη1ydy−1τd1∫x0x1∫0tΘ01x0−x,x1−yr1y,ςςm/τm+τmρdςdy−1τd1∑k=1∞b1,C1x1−x,sk1,τsk1/τd1Δ1′sk1∫x0x1∫0tesk1t−ς/τr1y,ςa1,C1x0−y,sk1,τςm/τm+τmρdςdy−1τd1∑k=1∞esk1t/τb1,C1x1−x,sk1,τsk1/τd1Δ1′sk1∫x0x1a1,C1x0−y,sk1,τη1ydy,(124)θ2x,t,τ=−1τd2∫x1x2Θ02x1−x,x2−yη2ydy−1τd2∫x1x2∫0tΘ02x1−x,x2−yr2y,ςςm/τm+τmρdςdy−1τd2∑k=1∞b2,C2x2−x,sk2,τsk2/τd2Δ2′sk2∫x1x2∫0tesk2t−ς/τr2y,ςa2,C2x1−y,sk2,τςm/τm+τmρdςdy−1τd2∑k=1∞esk2t/τb2,C2x2−x,sk2,τsk2/τd2Δ2′sk2∫x1x2a2,C2x1−y,sk2,τη2ydy,whereΘ01andΘ02 can be defined as in Lemma 4.Also, similar to (93), with the respective forms Φ01andΦ02 from Lemma 5 and the matching condition ξ1t=ζ2t, we get
(125)φ1x,t,τ=λt,τζ1tΦ01x1−x;b1+∑k=1∞Γk1ζ1tsk1Δ1′sk1b1,C1x1−x,sk1,τ−λt,τξ1tΦ01x0−x;a1−∑k=1∞Γk1ξ1tsk1Δ1′sk1a1,C1x0−x,sk1,τ+θ1x,t,τ,(126)φ2x,t,τ=λt,τξ1tΦ02x2−x;b2+∑k=1∞Γk2ξ1tsk2Δ2′sk2b2,C2x2−x,sk2,τ−λt,τξ2tΦ02x1−x;a2−∑k=1∞Γk2ξ2tsk2Δ2′sk2a2,C2x1−x,sk2,τ+θ2x,t,τ,where the operators Γk1andΓk2 are obtained from (94). The matching condition φ1x1,t=Λ1φ2x1,t yields
(127)Λ1λt,τΦ02x2−x1;b2+λt,τΦ01x0−x1;a1ξ1t+∑k=1∞a1,C1x0−x1,sk1,τsk1Δ1′sk1Γk1ξ1t−Λ1b2,C2x2−x1,sk2,τsk2Δ2′sk2Γk2ξ1t=λt,τζ1tΦ010;b1+θ1x1,t,τ+∑k=1∞b1,C10,sk1,τsk1Δ1′sk1Γk1ζ1t+Λ1λt,τξ2tΦ020;a2+∑k=1∞Λ1a2,C20,sk2,τsk2Δ2′sk2Γk2ξ2t−Λ1θ2x1,t,τ.For the unknown functionξ1, we can rewrite the linear integral equation (127) as
(128)λt,τξ1t+∑k=1∞akΓk1ξ1t+bkΓk2ξ1t=ct,where
(129)ak=a1,C1x0−x1,sk1,τsk1Δ1′sk1Λ1Φ02x2−x1;b2+Φ01x0−x1;a1,bk=−Λ1b2,C2x2−x1,sk2,τsk2Δ2′sk2Λ1Φ02x2−x1;b2+Φ01x0−x1;a1,ct=1Λ1Φ02x2−x1;b2+Φ01x0−x1;a1λt,τζ1tΦ010;b1+θ1x1,t,τ+∑k=1∞b1,C10,sk1,τsk1Δ1′sk1Γk1ζ1t+Λ1λt,τξ2tΦ020;a2+∑k=1∞Λ1a2,C20,sk2,τsk2Δ2′sk2Γk2ξ2t−Λ1θ2x1,t,τ.Inspired by the convolution formula (31), the natural transform of (128) is
(130)ℕξ1tλt,τs,τ+∑k=1∞akℕξ1tλt,τs,τ−τsk1s−sk1ℕξ1tλt,τ+bkℕξ1tλt,τs,τ−τsk2s−sk2ℕξ1tλt,τ=ℕcts,τ,which can be rewritten as
(131)1−ℕψts,τℕξ1tλt,τs,τ=ℕcts,τ.That is,(132)ℕξ1tλt,τs,τ=11−ℕψts,τℕcts,τ=1+ℕψts,τ+ℕψts,τ2+⋯ℕcts,τ=ℕcts,τ+ℕcts,τ∑m=1∞ℕψts,τm,where
(133)ℕψts,τ=−∑k=1∞ak1−τsk1s−sk1+bk1−τsk2s−sk2,for which the inverse natural transform is
(134)ψt=−τ∑k=1∞ak+bkδ0t+τ∑k=1∞aksk1esk1/τt+bksk2esk2/τt,where δ0 is the Dirac delta function. Hence, we have
(135)ξ1t=1λt,τct+∑m=1∞ψm∗ct,where for m≥2,ψm is the m-times self-convolution of ψ. Thus, one can conclude the solution to the bilayer diffusion problem (112)–(119) by the formulas (125) and (126), together with (123) and (124), with ξ1 given in (135). Now, it is time to attack an illustrative example in the following subsection.


### 3.1.1. Illustrative Example

Here, we discuss the solvability of the following two-layer diffusion system.Example 3.
Temperature distribution in the two-layer slab with mixed boundary condition.
Consider the following initial boundary value system. The diffusion equations are(136)∂φk∂t=dk∂2φk∂x2,xk−1<x<xk,x0=0,t>0,k=1,2,the initial conditions are
(137)φkx,0=ηkx,xk−1≤x≤xk,k=1,2,the outer boundary conditions are
(138)φ10,t=0,∂φ2x2,t∂x=0,t≥0,and the interface conditions at x=x1 are
(139)φ1x1,t=φ2x1,t,d1∂φ1x1,t∂x=d2∂φ2x1,t∂x,t≥0.
Comparing this problem with the general one that is defined by (112)–(119) reveals
(140)a1=1,0,b1=0,d1,a2=0,d2,b2=0,1,r1x,t=r2x,t=ζ1t=ξ2t=0,ξ1t=ζ2t.Here, we will consider ρ=0andτ=1.
Then, the solution of the problem (136)–(139) can be obtained from (125) and (126) as
(141)φ1x,t=−ξ1tΦ01x0−x;a1−∑k=1∞Γk1ξ1tsk1Δ1′sk1a1,C1x0−x,sk1,τ+θ1x,t,(142)φ2x,t=ξ1tΦ02x2−x;b2+∑k=1∞Γk2ξ1tsk2Δ2′sk2b2,C2x2−x,sk2,τ+θ2x,t,with
(143)θ1x,t=−1d1∫x0x1Θ01x0−x,x1−yη1ydy−1d1∑k=1∞esk1tb1,C1x1−x,sk1sk1/d1Δ1′sk1∫x0x1a1,C1x0−y,sk1η1ydy,θ2x,t=−1d2∫x1x2Θ02x1−x,x2−yη2ydy−1d2∑k=1∞esk2tb2,C2x2−x,sk2sk2/d2Δ2′sk2∫x1x2a2,C2x1−y,sk2η2ydy.The operatorsΓk1andΓk2 are obtained from (94).Next, we are going to simplify these formulas. Direct computations give(144)Δ1s=d1sd1coshsd1x1,Δ2s=−d2sd2sinhsd2x2−x1,Δ1′s=d1x12d1sinhsd1x1+12d1scoshsd1x1,Δ2′s=−d21d2sinhsd2x2−x1+x2−x12d2sd2coshsd2x2−x1.It is clear thatΔ1sk1=0andΔ2sk2=0 when
(145)sk1=2k−12π2d14x12,sk2=−k2π2d2x2−x12,k=1,2,3,⋯,respectively. Moreover,
(146)a1,C1x;sb1,C1y;ss/d1Δ1s=sinhs/d1xcoshs/d1ys/d1coshs/d1x1=O1,a2,C2x;sb2,C2y;ss/d2Δ2s=−coshs/d2xcoshs/d2ys/d2sinhs/d2x2−x1=O1s,as s⟶0. Therefore, (70) gives Θ01x,y=0andΘ02x,y=−d2/x2−x1.Further,(147)a1,C1−y;sk1=−isin2k−1π2x1y,b1,C1x1−x;sk1=−1k+12k−1πd12x1isin2k−1π2x1x,sk1d1Δ1′sk1=−1k2k−1π4,k=1,2,3,⋯,a2,C2x1−y;sk2=kπd2ix2−x1coskπx2−x1x1−y,b2,C2x2−x;sk2=kπix2−x1coskπx2−x1x2−x,sk2d2Δ2′sk2=−1kk2π22x2−x1.Thus,(148)θ1x,t=∑k=1∞2x1∫0x1η1ysin2k−1πy2x1dye−2k−12π2/4x12d1tsin2k−1π2x1x,θ2x,t=−1x2−x1∫x1x2η2ydy+∑k=1∞−1k1x2−x1∫x1x2η2ycoskπx2−x1x1−ydye−k2π22π2/x2−x12d2tcoskπx2−x1x2−x1.Similarly, we have(149)a1,C1x0−x,ssΔ1s=sinhs/d1x0−xsd1s/d1coshs/d1x1=O1s,b2,C2x2−x,ssΔ2s=−s/d2coshs/d2x2−xsd2s/d1sinhs/d2x2−x1=O1s2.So, (88) gives
(150)Φ01x0−x;a1=−xd1,Φ02x2−x;b2=16d2x2−x13x2−x2−x2−x12.Hence, loading the quantities in (141) and (142) gives the solution to problem (136)–(139) as
(151)φ1x,t=xd1ξ1t+8x1π2d1∑k=1∞−1kΓk1ξ1t2k−12sin2k−1π2x1x+θ1x,t,φ2x,t=−3x2−x2−x2−x126d2x2−x1ξ1t+2x2−xd2π2∑k=1∞−1kΓk2ξ1tk2coskπx2−x1x2−x+θ2x,t.Applying the matching conditionφ1x1,t=φ2x1,t gives
(152)ξ1t+∑k=1∞akΓk1ξ1t+bkΓk2ξ1t=ct,where
(153)ak=−24d2x12k−12π23d2−d1x1+d1x2,bk=−6d1x2−x1k2π23d2−d1x1+d1x2,ct=3d1d23d2−d1x1+d1x2θ1x1,t−θ2x1,t.To solve this integral equation, we use (134), to obtain
(154)ψt=δ0t−6d1d23d2−d1x1+d1x2∑k=1∞1x1e−2k−12π2d1/4x12t+1x2−x1e−k2π2d2/x2−x12t,where δ0 is the Dirac delta function. Here, we used
(155)∑k=1∞12k−12=π28,∑k=1∞1k2=π26.Hence, from (135), we have
(156)ξ1t=ct+ψ∗ct+ψ∗ψ∗ct+⋯=ct+∑m=1∞ψm∗ct.


## 3.1.1. Illustrative Example

Here, we discuss the solvability of the following two-layer diffusion system.Example 3.
Temperature distribution in the two-layer slab with mixed boundary condition.
Consider the following initial boundary value system. The diffusion equations are(136)∂φk∂t=dk∂2φk∂x2,xk−1<x<xk,x0=0,t>0,k=1,2,the initial conditions are
(137)φkx,0=ηkx,xk−1≤x≤xk,k=1,2,the outer boundary conditions are
(138)φ10,t=0,∂φ2x2,t∂x=0,t≥0,and the interface conditions at x=x1 are
(139)φ1x1,t=φ2x1,t,d1∂φ1x1,t∂x=d2∂φ2x1,t∂x,t≥0.
Comparing this problem with the general one that is defined by (112)–(119) reveals
(140)a1=1,0,b1=0,d1,a2=0,d2,b2=0,1,r1x,t=r2x,t=ζ1t=ξ2t=0,ξ1t=ζ2t.Here, we will consider ρ=0andτ=1.
Then, the solution of the problem (136)–(139) can be obtained from (125) and (126) as
(141)φ1x,t=−ξ1tΦ01x0−x;a1−∑k=1∞Γk1ξ1tsk1Δ1′sk1a1,C1x0−x,sk1,τ+θ1x,t,(142)φ2x,t=ξ1tΦ02x2−x;b2+∑k=1∞Γk2ξ1tsk2Δ2′sk2b2,C2x2−x,sk2,τ+θ2x,t,with
(143)θ1x,t=−1d1∫x0x1Θ01x0−x,x1−yη1ydy−1d1∑k=1∞esk1tb1,C1x1−x,sk1sk1/d1Δ1′sk1∫x0x1a1,C1x0−y,sk1η1ydy,θ2x,t=−1d2∫x1x2Θ02x1−x,x2−yη2ydy−1d2∑k=1∞esk2tb2,C2x2−x,sk2sk2/d2Δ2′sk2∫x1x2a2,C2x1−y,sk2η2ydy.The operatorsΓk1andΓk2 are obtained from (94).Next, we are going to simplify these formulas. Direct computations give(144)Δ1s=d1sd1coshsd1x1,Δ2s=−d2sd2sinhsd2x2−x1,Δ1′s=d1x12d1sinhsd1x1+12d1scoshsd1x1,Δ2′s=−d21d2sinhsd2x2−x1+x2−x12d2sd2coshsd2x2−x1.It is clear thatΔ1sk1=0andΔ2sk2=0 when
(145)sk1=2k−12π2d14x12,sk2=−k2π2d2x2−x12,k=1,2,3,⋯,respectively. Moreover,
(146)a1,C1x;sb1,C1y;ss/d1Δ1s=sinhs/d1xcoshs/d1ys/d1coshs/d1x1=O1,a2,C2x;sb2,C2y;ss/d2Δ2s=−coshs/d2xcoshs/d2ys/d2sinhs/d2x2−x1=O1s,as s⟶0. Therefore, (70) gives Θ01x,y=0andΘ02x,y=−d2/x2−x1.Further,(147)a1,C1−y;sk1=−isin2k−1π2x1y,b1,C1x1−x;sk1=−1k+12k−1πd12x1isin2k−1π2x1x,sk1d1Δ1′sk1=−1k2k−1π4,k=1,2,3,⋯,a2,C2x1−y;sk2=kπd2ix2−x1coskπx2−x1x1−y,b2,C2x2−x;sk2=kπix2−x1coskπx2−x1x2−x,sk2d2Δ2′sk2=−1kk2π22x2−x1.Thus,(148)θ1x,t=∑k=1∞2x1∫0x1η1ysin2k−1πy2x1dye−2k−12π2/4x12d1tsin2k−1π2x1x,θ2x,t=−1x2−x1∫x1x2η2ydy+∑k=1∞−1k1x2−x1∫x1x2η2ycoskπx2−x1x1−ydye−k2π22π2/x2−x12d2tcoskπx2−x1x2−x1.Similarly, we have(149)a1,C1x0−x,ssΔ1s=sinhs/d1x0−xsd1s/d1coshs/d1x1=O1s,b2,C2x2−x,ssΔ2s=−s/d2coshs/d2x2−xsd2s/d1sinhs/d2x2−x1=O1s2.So, (88) gives
(150)Φ01x0−x;a1=−xd1,Φ02x2−x;b2=16d2x2−x13x2−x2−x2−x12.Hence, loading the quantities in (141) and (142) gives the solution to problem (136)–(139) as
(151)φ1x,t=xd1ξ1t+8x1π2d1∑k=1∞−1kΓk1ξ1t2k−12sin2k−1π2x1x+θ1x,t,φ2x,t=−3x2−x2−x2−x126d2x2−x1ξ1t+2x2−xd2π2∑k=1∞−1kΓk2ξ1tk2coskπx2−x1x2−x+θ2x,t.Applying the matching conditionφ1x1,t=φ2x1,t gives
(152)ξ1t+∑k=1∞akΓk1ξ1t+bkΓk2ξ1t=ct,where
(153)ak=−24d2x12k−12π23d2−d1x1+d1x2,bk=−6d1x2−x1k2π23d2−d1x1+d1x2,ct=3d1d23d2−d1x1+d1x2θ1x1,t−θ2x1,t.To solve this integral equation, we use (134), to obtain
(154)ψt=δ0t−6d1d23d2−d1x1+d1x2∑k=1∞1x1e−2k−12π2d1/4x12t+1x2−x1e−k2π2d2/x2−x12t,where δ0 is the Dirac delta function. Here, we used
(155)∑k=1∞12k−12=π28,∑k=1∞1k2=π26.Hence, from (135), we have
(156)ξ1t=ct+ψ∗ct+ψ∗ψ∗ct+⋯=ct+∑m=1∞ψm∗ct.


## 3.2. Solution of a Multilayer Problem

Here, we investigate the solvability of the main problem (2)–(8), through solving the initial boundary value problems (9)–(11). Similar to what we have denoted in Section 2, we consider the following notations:
(157)aj=ı,ι,j=1,vj,μj,2≤j≤n,bj=vj,μj,1≤j≤n−1,ℓ,l,j=n,and, for all 1≤j≤n,
(158)Ljy;s,τ=coshsτdjy,sτdjsinhsτdjy,Cjy;s,τ=sinhsτdjy,sτdjcoshsτdjy.Moreover, we define(159)Δjs=bj,Cjxj,s,τaj,Ljxj−1,s,τ−aj,Cjxj−1,s,τbj,Ljxj,s,τ,and let skjk=1∞ be the sequence of zeros of the function Δjsforall1≤j≤n, i.e.,
(160)Δjskj=0,Δj′skj≠0k=1,2,⋯.Analogue to the computations of (84) and (93), we have for the current case, for all j=1,⋯,n,
(161)θjx,t,τ=−1τdj∫xj−1xjΘ0jxj−1−x,xj−yηjydy−1τdj∫xj−1xj∫0tΘ0jxj−1−x,xj−yrjy,ζζm/τm+τmρdζdy−1τdj∑k=1∞bj,Cjxj−x,skj,τskj/τdjΔj′skj∫xj−1xj∫0teskjt−ζ/τrjy,ζaj,Cjxj−1−y,skj,τζm/τm+τmρdζdy−1τdj∑k=1∞eskjt/τbj,Cjxj−x,skj,τskj/τdjΔj′skj∫xj−1xjaj,Cjxj−1−y,skj,τηjydy,where Θ0j can be defined in a similar way as in Lemma 4, and
(162)φjx,t,τ=λt,τζjtΦjxj−x,0;bj−∫0tλς,τD2Φjxj−x,t−ς;bjζjςdς−λt,τξjtΦjxj−1−x,0;aj+∫0tλς,τD2Φjxj−1−x,t−ς;ajξjςdς+θjx,t,τ,with the respective forms Φj defined by (87) in Lemma 5. This last equation (162) can be rewritten as
(163)φjx,t,τ=Tjζ~jxj−x,t;bj−Tjξ~jxj−1−x,t;aj+θjx,t,τ,in which ζ~j=λt,τζjt,ξ~j=λt,τξjt and the linear operator Tj is defined by
(164)Tjφy,t;L=Φjy,0;Lφt−∫0tD2Φjy,t−ς;Lφςdς,for all j=1,⋯,n,L∈ℝ2.The matching conditionsφjxj,t,τ=Λjφj+1xj,t,τ,j=1,⋯,n−1, lead to
(165)Tjζ~j0,t;bj−Tjξ~jxj−1−xj,t;aj+θjxj,t,τ=ΛjTj+1ζ~j+1xj+1−xj,t;bj+1−Tj+1ξ~j+10,t;aj+1+θj+1xj,t,τ.Using the matching conditions (12), we have ζ~j+1t=ξ~jt for all 1≤j≤n−1. Thus, for j=1,
(166)−T1ξ~1x0−x1,t;a1−Λ1T2ξ~1x2−x1,t;b2+Λ1T2ξ~20,t;a2=Λ1θ2x1,t,τ−θ1x1,t,τ−T1ζ~10,t;b1.For2≤j≤n−2,
(167)Tjξ~j−10,t;bj−Tjξ~jxj‐1−xj,t;aj−ΛjTj+1ξ~jxj+1−xj,t;bj+1+ΛjTj+1ξ~j+10,t;aj+1=Λjθj+1xj,t,τ−θjxj,t,τ.Forj=n−1,
(168)Tn−1ξ~n−20,t;bn−1−Tn−1ξ~n−1xn−2−xn−1,t;an−1−Λn−1Tnξ~n−1xn−xn−1,t;bn=Λn−1θnxn−1,t,τ−θn−1xn−1,t,τ−Λn−1Tnξ~n0,t;an.System (166), (167), and (168), of n−1 integral equations of the unknowns ξ~j;1≤j≤n−1, can be adjusted as a matrix equation
(169)A0ht+A′∗ht=bt,with At is a tridiagonal matrix of order n−1 whose entries are as follows:
(170)−Φjxj−1−xj,t;aj−ΛjΦjxj+1−xj,t;bj+1,1⩽j⩽n−1maindiagonal,ΛjΦj+10,t;aj+1,1⩽j⩽n−2superdiagonalΦj0,t;bj,2⩽j⩽n−1subdiagonal,and the vectors ht and bt are defined as
(171)ht=ξ~1t⋮ξ~n−1t,bt=Λ1θ2x1,t,τ−θ1x1,t,τ−T1ζ~10,t;b1Λ2θ3x2,t,τ−θ2x2,t,τ⋮Λn−2θn−1xn−2,t,τ−θn−2xn−2,t,τΛn−1θnxn−1,t,τ−θn−1xn−1,t,τ−Λn−1Tnξ~n0,t;an.In fact, we can rewrite (169) as
(172)ht=Ct+B∗ht,with Ct=A0−1bt and Bt=−A0−1A′t. In view of the convolution formula (31), the natural transform of (172) reads
(173)ℕhts,τ=ℕCts,τ+τℕBts,τℕhts,τ,which is equivalent to
(174)ℕhts,τ=I−τℕBts,τ−1ℕCts,τ=I+τℕBts,τ+τℕBts,τ2+⋯ℕCts,τ,where I is the n−1×n−1 identity matrix. Once again, throughout the convolution sense (31), the natural transform inversion of (174) is
(175)ht=Ct+B∗Ct+B∗B∗Ct+⋯=Ct+∑m=1∞Bm∗Ct,where Bm is the m-times self-convolution of B. Finally, the solution of the nonhomogeneous multilayer diffusion systems (9)–(11) and hence that of the main problem (2)–(8) is concluded as
(176)φjx,t,τ=λt,τζjtΦ0jxj−x;bj−λt,τξjtΦ0jxj−1−x;aj+∑k=1∞λt,τζjt−∫0tλς,τskjeskjt−ς/τζjςdςskj∆j′skjbj,Cjxj−x,skj,τ−∑k=1∞λt,τξjt−∫0tλς,τskjeskjt−ς/τξjςdςskj∆j′skjaj,Cjxj−1−x,skj,τ+θjx,t,τ,with the respective forms Φ0j and θj defined as in (88) and (161), respectively, for all j=1,⋯,n.


## 4. Conclusion

Throughout the current contribution, a one-dimensionaln-layer nonhomogeneous diffusion problem with time-varying data and general interface conditions has been concluded by means of a generalized integral transform. Although most of the previous works have been focused on solving the problems of the homogeneous diffusion equation, the nonhomogeneous diffusion equation problem arises in many physical applications. We have obtained the exact solutions for one- and multilayer nonhomogeneous diffusion problems. The former case has been solved by a new generalized integral transform; the later one (n-layer problem) has been recast in a sequence of one-layer problems. The obtained results generalize and extend those in [31, 33, 38, 40, 42]. Our results motivate to deal with other types of diffusion problems, for example, reaction diffusion problems, advection-reaction diffusion problems, and nonautonomous reaction diffusion problems.On the other hand, more general partial differential equations (PDEs) and systems can be considered, for example, system of coupled PDEs, nonlinear diffusion PDEs, and nonautonomous reaction diffusion PDEs. Those kinds of PDEs appear widely as epidemiological models to study and analyze the spread of diseases and pandemics [22–25].



---

*Source: 2304219-2022-03-07.xml*