yo=[6*10^4 10^6 0 10^4 0 4.07 3.33 25.75 1.45 0 0 36.2 60 0 781 105 100 1.25 13.5 20 7.4 82.5];
MAPdata=[115 120 117 112 112.5 113 113 113 113 113 113 113 113];
if (tspan<=max(tsimulation))
MAP(:,k)=interp1(tsimulation,MAPdata,tspan,"linear");
MAP(:,k)=interp1(tsimulation,MAPdata,tspan,"linear","extrap");
[t(:,k),y(:,:,k),MAP,Nstep,xNOL]=ode45(@(t,y,MAP,Nstep,xNOL)ActualmodelwithPinfinity(t,y,MAP(:,k),Nstep,xNOL(k)),tspan,yo);
function dydt=ActualmodelwithPinfinity(t,y,MAP,Nstep,xNOL)
dydt(1)=kpg*y(1)-kd*y(1)-((sm*y(1))/(um+kpm*y(1)))-a*HU2(y(1),kcNA,1)*y(3)-b*HU2(y(1),kcMA,1)*y(5);
dydt(2)=knr*y(2)*(1-(y(2)/Ns))-kNP*HU2(y(1),xNP,hNP)*y(2)*(1+kNTNF*HU2(y(6),xNTNF,2))*(1+kN6*HU2(y(7),xN6,2))*(1+kN8*HU2(y(8),xN8,3))*HYX(y(9),xN10,2)-kNRd*y(2);
dydt(3)=kNP*HU2(y(1),xNP,hNP)*y(2)*(1+kNTNF*HU2(y(6),xNTNF,2))*(1+kN6*HU2(y(7),xN6,2))*(1+kN8*HU2(y(8),xN8,3))*HYX(y(9),xN10,2)-kNAd*y(3);
dydt(4)=kMR*(1-(y(4)/Minfinty))*y(4)-kMP*HU2(y(1),xMP,hMP)*y(4)*(1+kMTNF*HU2(y(6),xMTNF,2))*(1+kM6*HU2(y(7),xM6,2))*HYX(y(9),xM10,2)-kMRd*y(4);
dydt(5)=kMP*HU2(y(1),xMP,hMP)*y(4)*(1+kMTNF*HU2(y(6),xMTNF,2))*(1+kM6*HU2(y(7),xM6,2))*HYX(y(9),xM10,2)-kMAd*y(5);
dydt(6)=(kTNFN*HU2(y(3),xTNFN,hTNFNM)+kTNFM*HU2(y(5),xTNFM,hTNFNM))*HU2(y(16),xGTNF,2)*HYX(y(18),xLTNF,1)*HYX(y(7),xTNF6,2)*HYX(y(9),xTNF10,2)-kTNF*(y(6)-wTNF);
dydt(7)=(k6N*HU2(y(3),x6N,h6NM)+k6M*HU2(y(5),x6M,h6NM))*(k6TNF*HU2(y(6),x6TNF,2)+k6NO*HU2(y(11),x6N0,2))*HYX(y(18),xL6,1)*HU2(y(16),xG6,2)*HYX(y(7),x66,1)*HYX(y(9),x610,2)-k6*(y(7)-wIL6);
dydt(8)=(k8N*HU2(y(3),x8N,h8N)+k8TNF*HU2(y(6)-wTNF,x8TNF,3))*HYX(y(9),x810,2)-k8*(y(8)-wIL8);
dydt(9)=(k10N*HU2(y(3),x10N,4)+k10M*HU2(y(5),x10M,4))*(k10TNF*HU2(y(6),x10TNF,2)+k106*HU2(y(7),x106,4))-k10*(y(9)-wIL10);
dydt(10)=(u/(1+kA*y(9)))+kdp*hV(y(1),xdp)-ud*y(10);
dydt(11)=(kNONM*(y(3)+y(5))+kNOD*y(10))*HU2(y(6)-wTNF,nNOTNF,2)*HYX(y(9)-wIL10,nNO10,1)-uno*y(11);
dydt(12)=kTTNF*HU2((y(6)-wTNF),nTTNF,hT)+kT6*HU2((y(7)-wIL6),nT6,hT)-kT10*(1-HYX((y(9)-wIL10),nT10,hT))-kT*(y(12)-Tb);
fBP(j)=HU2(BPb-MAP(:,j),nHPlBP,hHP);
fBP(j)=HYX(BPb-MAP(:,j),nHPhBP,hHP);
dydt(13)=(Hb-y(13)+(kH*(HM-Hb)*HU2((y(12)-Tb),nHT,2)*fBP(j)))/tauH;
dydt(14)=kpe*HU2(y(1),npe,hET)-kde*y(14);
dydt(15)=-kPTe*y(14)*y(15)-kPT*(y(15)-PTb);
dydt(16)=(sg+kLG*(y(18)-Lact0)+(kglyglu*y(17)*(1+HU2(y(6)-wTNF,xTNFG,2))*(1+HU2(y(7)-wIL6,xIL6G,2))))*HYX(y(19)-wINSULIN,xING,1)-kglucgly*y(16)-sIN*y(16)*y(19);
dydt(17)=kglucgly*y(16)-kglyglu*y(17)*(1+HU2(y(6)-wTNF,xTNFG,2))*(1+HU2(y(7)-wIL6,xIL6G,2))*HYX(y(19)-wINSULIN,xING,1);
dydt(18)=kGL*y(16)*(1-HYX(y(11),xNOL,2))*(1+HU2(y(14),xETG,1))-uL*(y(18)-Lact0);
dydt(19)=kING*(1+HU2(y(16)-Go,alpha,2))*(1-HYX(y(6)-wTNF,xINTNF,1))*(1-HYX(y(7)-wIL6,xIN6,1))-kIN*(y(19)-wINSULIN);
dydt(20)=(RRb-y(20)+(kRRlact*(1+HU2(y(18),nRRlact,1))*HU2(y(12)-Tb,nTempRR,2)*HYX(y(22),nPaO2RR,1)))/tauRR;
dydt(21)=-kPHlact*y(18)*y(21)*(1-HYX(Pao2b-y(22),nPHPao2,1))+kpH*(PHbaseline-y(21));
dydt(22)=-kinfPao2*HU2(y(12)-Tb,nTempPao2,2)*HYX(y(21),nPaO2Ph,1)-ko2Pmetab*(y(22)-Pao2b);