Solving the steady-state problem of a multilayer hollow sphere

Question

I would like to reproduce the analytical solution for the steady state problem as described in the paper Singh2016(https://doi.org/10.1115/1.4033536).

The additional conditions can be found in the following implementation attempt, where I try to reproduce a T(r) result plot from the paper, for the time being only for the inner layer (i=1) and for the first 2 steps of the generalized Fourier series.

ClearAll[all];
r0=1;qheat0=1;k1=1;\[Alpha]1=1;k2=2 k1;k3=4 k1;\[Alpha]2=4 \[Alpha]1;\[Alpha]3=9 \[Alpha]1;r1=2 r0;r2=4 r0;r3=6r0;Ain=0;Bin=1;Aout=k3;Bout=0;Cin=0; \[Theta]=0;\[Phi]=0;
T1[r_,\[Theta]_,\[Phi]_]:=Sum[Tc1mp[r,\[Theta],\[Phi],m,q]*LegendreP[m,q,Cos[\[Theta]]]*Cos[q \[Phi]],{m,0,2},{q,0,m}]+Sum[Ts1mp[r,\[Theta],\[Phi],m,q]*LegendreP[m,q,Cos[\[Theta]]]*Sin[q \[Phi]],{m,1,2},{q,1,m}];

Tc1mp[r_,\[Theta]_,\[Phi]_,m_,q_]:=Ac1mp[r,\[Theta],\[Phi],m,q] r^m+Bc1mp[r,\[Theta],\[Phi],m,q] r^(-m-1);\[IndentingNewLine]Ts1mp[r_,\[Theta]_,\[Phi]_,m_,q_]:=As1mp[r,\[Theta],\[Phi],m,q] r^m+Bs1mp[r,\[Theta],\[Phi],m,q] r^(-m-1);
\[IndentingNewLine]{{As1mp[r_,\[Theta]_,\[Phi]_,m_,q_]},{Bs1mp[r_,\[Theta]_,\[Phi]_,m_,q_]},{As2mp[r_,\[Theta]_,\[Phi]_,m_,q_]},{Bs2mp[r_,\[Theta]_,\[Phi]_,m_,q_]},{As3mp[r_,\[Theta]_,\[Phi]_,m_,q_]},{Bs3mp[r_,\[Theta]_,\[Phi]_,m_,q_]}}:=Inverse[{{(Ain m+Bin r0) r0^(m-1),(-Ain(m+1)+Bin r0) r0^(-m-2),0,0,0,0},{r1^m,r1^(-m-1),-r1^m,-r1^(-m-1),0,0},{0,0,r2^m,r2^(-m-1),-r2^m,-r2^(-m-1)},{k1 m r1^(m-1),-k1 (m+1) r1^(-m-2),-k2 m r1^(m-1),k2(m+1) r1^(-m-2),0,0},{0,0,k2 m r2^(m-1),(-k2)(m+1)r2^(-m-2),-k3 m r2^(m-1),k3(m+1) r2^(-m-2)},{0,0,0,0,(Aout m+Bout r2) r3^(m-1),(-Aout(m+1)+Bout r3) r3^(-m-2)}}]. {{0},{0},{0},{0},{0},{Coutsmp[\[Theta],\[Phi],m, q]}};

{{Ac1mp[r_,\[Theta]_,\[Phi]_,m_,q_]},{Bc1mp[r_,\[Theta]_,\[Phi]_,m_,q_]},{Ac2mp[r_,\[Theta]_,\[Phi]_,m_,q_]},{Bc2mp[r_,\[Theta]_,\[Phi]_,m_,q_]},{Ac3mp[r_,\[Theta]_,\[Phi]_,m_,q_]},{Bc3mp[r_,\[Theta]_,\[Phi]_,m_,q_]}}:=Inverse[{{(Ain m+Bin r0) r0^(m-1),(-Ain(m+1)+Bin r0) r0^(-m-2),0,0,0,0},{r1^m,r1^(-m-1),-r1^m,-r1^(-m-1),0,0},{0,0,r2^m,r2^(-m-1),-r2^m,-r2^(-m-1)},{k1 m r1^(m-1),-k1 (m+1) r1^(-m-2),-k2 m r1^(m-1),k2(m+1) r1^(-m-2),0,0},{0,0,k2 m r2^(m-1),(-k2)(m+1)r2^(-m-2),-k3 m r2^(m-1),k3(m+1) r2^(-m-2)},{0,0,0,0,(Aout m+Bout r2) r3^(m-1),(-Aout(m+1)+Bout r3) r3^(-m-2)}}]. {{0},{0},{0},{0},{0},{Coutcmp[\[Theta],\[Phi], m, q]}};

Zq[q_]:=Piecewise[{{\[Pi],q!=0},{2\[Pi],q=0}}];
Nmq[m_,q_]:=(2/(2m+1)) (m+q)!/(m-q)!;

Coutcmp[\[Theta]_,\[Phi]_,m_,q_]:=(1/(Nmq[m,q ] Zq[q])) Integrate[qheat[\[Theta],\[Phi]]*LegendreP[m,q,Cos[\[Theta]]] *Cos[q \[Phi]],{\[Phi],0,2\[Pi]},{\[Theta],0,\[Pi]}];
Coutsmp[\[Theta]_,\[Phi]_,m_,q_]:=(1/(Nmq[m,q] \[Pi]))   Integrate[qheat[\[Theta],\[Phi]]*LegendreP[m,q,Cos[\[Theta]]]* Sin[q \[Phi]],{\[Phi],0,2\[Pi]},{\[Theta],0,\[Pi]}];

qheat[\[Theta]_,\[Phi]_]:=If[\[Phi]<\[Pi],(qheat0 \[Theta]^2 (\[Pi]-\[Theta])^2 \[Phi]^2 (\[Pi]-\[Phi])^2),0];
\[IndentingNewLine]sol = SolveValue[{T1[r,\[Theta],\[Phi]]},r]
Plot[sol,{r,1,2}]

Unfortunately I have difficulties to initialize the functions which are set up by the matrix I would be grateful for any advice!

Edit: Entirely new set-up