I’m currently trying to replicate the process covered in a paper dealing with a transient heat transfer setup in a multi-layered cylinder. (https://doi.org/10.1016/j.ijheatmasstransfer.2016.12.020)

I’ve gotten a good chunk of my code to work properly as I have gotten the same roots/eigenvalues as the paper got for the 2-layer cylinder example. My issue lies in the portion after this as the only thing I can compare to the paper is my final result.

Below is my code, the expected graph, and the resulting graph from my code in that order

import scipy as sci
import numpy as np
import matplotlib.pyplot as plt
from chebpy import chebfun

'''
inputs!
radius layers
thermal conductivities (k)
thermal diffusivities (alpha)
inner surface temp
outer surface temp
'''

radii = np.array([1,2,4])
k=np.array([4,1])
alpha=np.array([4,1])
Tin=1
Tout=0
Mb=np.array([[1,0],[0,0]])
Nb=np.array([[0,0],[1,0]])
'''
-----------------------------------------------------------------------------------------------------------------------
Eigenvalues and roots
'''
def lam(s,i): #between eq 21 and 22
    asqr=-s/alpha[i]
    a=asqr**0.5
    return a

def i_check(r): #gives an i value based on radius
    ival=1
    check=0
    while check==0:
        if r<=radii[ival]:
            check=1
        else:
            ival+=1
    return ival

def u_anyr(r,s,tog): #gives u or inv u for any r and s, if tog is 0, gives u, otherwise inv u, eq 23
    i=i_check(r)
    if tog==0:
        u = np.array([[sci.special.j0(lam(s, i - 1) * r), sci.special.y0(lam(s, i - 1) * r)],
                      [-1 * lam(s, i - 1) * k[i - 1] * sci.special.j1(lam(s, i - 1) * r), -1 * lam(s, i - 1) * k[i - 1] * sci.special.y1(lam(s, i - 1) * r)]])
    else: # auto jumps to previous layer
        u = np.array([[sci.special.j0(lam(s, i - 1) * radii[i-1]), sci.special.y0(lam(s, i - 1) * radii[i-1])],
                      [-1 * lam(s, i - 1) * k[i - 1] * sci.special.j1(lam(s, i - 1) * radii[i-1]),
                       -1 * lam(s, i - 1) * k[i - 1] * sci.special.y1(lam(s, i - 1) * radii[i-1])]])
        u=np.linalg.inv(u)
    return u
def mu_anyr(r,s):#gives mu for any r, eq 30
    return np.matmul(u_anyr(r,s,0),u_anyr(r,s,1))
def U_anyr(r,s): #gives U for any r, eq 29
    i=i_check(r)
    q=np.array([[1,0],[0,1]])
    check = 0
    layer=i
    while check==0:
        if layer == i:
            q=np.matmul(q,mu_anyr(r,s))
        else:
            q=np.matmul(q,mu_anyr(radii[layer],s))
        if layer==1:
            check=1
        layer-=1
    return q

def Z(s):
    U=U_anyr(radii[-1],s)
    n=np.matmul(Nb,U)
    m=Mb+n
    return m

def f(x):
    a=np.array([])
    for j in range(len(x)):
        a = np.append(a,np.linalg.det(Z(x[j])))
    return a
'''
x=np.linspace(-100,10,1000)
y=f(x)
plt.plot(x,y)
plt.show()
'''
f = chebfun(lambda x: f(x),[-1000000,-0.001])
eigen = f.roots()
negbsquare=np.flip(eigen)
print(negbsquare.size)
def de_neg_and_square(num):
    return (num*-1)**0.5
just_root=np.vectorize(de_neg_and_square)
bees=just_root(negbsquare)
r=negbsquare[0]
'''
-----------------------------------------------------------------------------------------------------------------------
Residues
'''
def d_mu(i,s):
    qw = (s*-1)**0.5
    du = np.array([[((-1*radii[i])/(alpha[i-1]**0.5))*sci.special.j1(lam(s,i-1)*radii[i]),
                    ((-1*radii[i])/(alpha[i-1]**0.5))*sci.special.y1(lam(s,i-1)*radii[i])],
                   [((-1*qw*k[i-1]*radii[i])/(alpha[i-1]))*sci.special.j0(lam(s,i-1)*radii[i]),
                    ((-1*qw*k[i-1]*radii[i])/(alpha[i-1]))*sci.special.y0(lam(s,i-1)*radii[i])]])

    st=(radii[i-1]*np.pi)/(2*k[i-1])
    inv_du = st* np.array([[((-1*qw*k[i-1]*radii[i-1])/(alpha[i-1]))*sci.special.y0(lam(s,i-1)*radii[i-1]),
                            ((radii[i-1])/(alpha[i-1]**0.5))*sci.special.y1(lam(s,i-1)*radii[i-1])],
                           [((qw*k[i-1]*radii[i-1])/(alpha[i-1]))*sci.special.j0(lam(s,i-1)*radii[i-1]),
                            ((-1*radii[i-1])/(alpha[i-1]**0.5))*sci.special.j1(lam(s,i-1)*radii[i-1])]])
    a=np.matmul(du,u_anyr(radii[i],s,1))
    b=np.matmul(u_anyr(radii[i],s,0),inv_du)
    p=a+b
    q=-1/(2*qw)
    qwer=q*p
    return qwer
#print(d_mu(1,r))
def dU(s):
    aw=np.array([[0.0,0.0],[0.0,0.0]])
    i=i_check(radii[-1])
    for j in range(len(alpha)):
        q=np.array([[1,0],[0,1]])
        for h in range(len(alpha)):
            if h == i-1:
                q=np.matmul(d_mu(i,s),q)
            else:
                q=np.matmul(mu_anyr(radii[h+1],s),q)
        aw+=q
        i-=1
    return aw
#print(dU(r))
def Resi(s):
    b=(s*-1)**0.5
    z=Z(s)
    a=np.array([[0.0,0.0],[0.0,0.0]])
    a[0,0]=z[1,1]
    a[1,1]=z[0,0]
    a[1,0]=-1*z[1,0]
    a[0,1]=-1*z[0,1]
    n=(1/(2*b))*np.matmul(Nb,dU(s))
    m=Mb-n
    d=np.linalg.det(m)
    return a/d


#print(np.matmul(Resi(r),Z(r)))
'''
------------------------------------------------------------------------------------------------------------------------
Answer?
'''
def Ibk(t,s):
    f1 = lambda x:np.exp(s*(t-x))*Tin
    f2 = lambda x:np.exp(s*(t-x))*Tout
    If1,err = sci.integrate.quad(f1,0,t)
    If2,err = sci.integrate.quad(f2,0,t)
    return np.array([[If1],[If2]])
def n(r,t,s):
    a=np.array([[0.0],[0.0]])
    op=np.zeros((len(r),2,1))
    for l in range(len(r)):
        for j in range(len(s)):
            qe = U_anyr(r[l], s[j])
            qe = np.matmul(qe, Resi(s[j]))
            qe = np.matmul(qe, Ibk(t, s[j]))
            a += qe
        op[l]=a
    return op


#print(n([2],0.1,negbsquare))

x=np.linspace(radii[0],radii[-1],30)
y=n(x,0.8,negbsquare)[:,0,:]
plt.plot(x,y)
plt.show()

The wanted graph
enter image description here

The resulting graph
enter image description here