Adaptive step Size ODE45

Hi! I am currently trying to implement the adaptive step size method to solve my system of ODEs. I have coded out my ODEs based on the RK 4th order method but however with a Fixed step size of 1 Second. I was hoping to adopt adaptive step sizes so as to yield more accurate results. Any help or advice at all will be greatly appreciated! Is there any way for me to loop the step sizes to get adaptive step sizes?

Here is my code below.

if true
  format long
%to insert upper boundary conditions here%
T01=37.404;
T02=36.235;
T03=36.017;
T04=35.806;
T05=37.750;
T06=36.839;
T07=33.499;
T08=31.834;
T09=36.039;
T010=34.850;
T011=33.666;
T012=33.162;
T013=33.330;
T014=33.142;
T015=33.024;
T016=32.859;
T017=33.446;
T018=32.932;
T019=31.895;
T020=31.666;
T021=31.528;
T022=31.477;
T023=31.428;
T024=31.345;
T025=37.092;

%end of up top boundary conditions

y1=T01;
y2=T02;
y3=T03;
y4=T04;
y5=T05;
y6=T06;
y7=T07;
y8=T08;
y9=T09;
y10=T010;
y11=T011;
y12=T012;
y13=T013;
y14=T014;
y15=T015;
y16=T016;
y17=T017;
y18=T018;
y19=T019;
y20=T020;
y21=T021;
y22=T022;
y23=T023;
y24=T024;
y25=T025;
%end of upper boundary allocation

a=0; 
b=8; %time interval, 0 to 8 hrs
Dt=(1)/3600; %Fixed Step Size
N=(b-a)/Dt; %steps

%t0=29.7; %operative temp, set to 29.7 from manikin experiment papers%

   %Preallocate arrays
  MT1=zeros(1,N);
  MT2=zeros(1,N);
  MT3=zeros(1,N);
  MT4=zeros(1,N);
  MT5=zeros(1,N);
  MT6=zeros(1,N);
  MT7=zeros(1,N);
  MT8=zeros(1,N);
  MT9=zeros(1,N);
  MT10=zeros(1,N);
  MT11=zeros(1,N);
  MT12=zeros(1,N);
  MT13=zeros(1,N);
  MT14=zeros(1,N);
  MT15=zeros(1,N);
  MT16=zeros(1,N);
  MT17=zeros(1,N);
  MT18=zeros(1,N);
  MT19=zeros(1,N);
  MT20=zeros(1,N);
  MT21=zeros(1,N);
  MT22=zeros(1,N);
  MT23=zeros(1,N);
  MT24=zeros(1,N);
  MT25=zeros(1,N);

t=zeros(1,N);
t(1)=a;
%end of preallocation

MT1(1)=T01;
MT2(1)=T02;
MT3(1)=T03;
MT4(1)=T04;
MT5(1)=T05;
MT6(1)=T06;
MT7(1)=T07;
MT8(1)=T08;
MT9(1)=T09;
MT10(1)=T010;
MT11(1)=T011;
MT12(1)=T012;
MT13(1)=T013;
MT14(1)=T014;
MT15(1)=T015;
MT16(1)=T016;
MT17(1)=T017;
MT18(1)=T018;
MT19(1)=T019;
MT20(1)=T020;
MT21(1)=T021;
MT22(1)=T022;
MT23(1)=T023;
MT24(1)=T024;
MT25(1)=T025;

Vmove=0;
Vair=0.2;
if Vmove == 0
    hcnew=11.6*(Vair^0.5);
  else
    hcnew=(8.6*(Vmove^0.53))+(1.96*(Vair*0.86));
end

Mhcnew=hcnew*area;

Tr=40; %mean radiant temperature assumed to be 40
ta=33; %ambient temp assumed to be 32 degrees
pa=1.9; %Vapor pressure at head segment, 1.9kpa?
%t0=((hcnew*ta)+(hr*Tr))./(hr+hcnew);
%t0=((hc*ta)+(hr*Tr))./(hr+hc);
t0=29.7; %operative temp, set to 29.7 from manikin experiment papers%

for i=1:(N-1)

T1=y1;
T2=y2;
T3=y3;
T4=y4;
T5=y5;
T6=y6;
T7=y7;
T8=y8;
T9=y9;
T10=y10;
T11=y11;
T12=y12;
T13=y13;
T14=y14;
T15=y15;
T16=y16;
T17=y17;
T18=y18;
T19=y19;
T20=y20;
T21=y21;
T22=y22;
T23=y23;
T24=y24;
T25=y25;

T=[T1;T2;T3;T4;T5;T6;T7;T8;T9;T10;T11;T12;T13;T14;T15;T16;T17;T18;T19;T20;T21;T22;T23;T24;T25];

Err=(T-Tset);

Wrm=max(Err,0);
Cld=min(Err,0)*(-1);
Wrmsold=SKINR.*Wrm(4:4:end,:); %replae i and j%
Wrms=sum(Wrmsold);

Cldsold=SKINR.*Cld(4:4:end,:);
Clds=sum(Cldsold);

ST=-(Cst)*Err(1)-Sst*(Wrms-Clds);

DL=Cdl*Err(1)+Sdl*(Wrms-Clds);

km=2.0.^(Err(4:4:end,:)/10);

Esw=(Csw*Err(1)+Ssw*(Wrms-Clds))*(SKINS.*km); 
%SKINS & km - 6x1 matrixes%
%Esw - 6x1 matrix%

E=Eb+Esw; %Evaporative heat loss at skin surface - 6x1 matrix%
%sum of evaporative heat loss of sweat and water vapor diffusion at skin layer%
%Values for Eb used from Stolwijk Model

Qt=(ht.*(T(4:4:end,:)-t0)).*area; 

%HEAD SEGMENT%
EK11= Dt * ((1/C(1,1))*(Qb(1,1)-pC*BFB(1,1)*(T1-T25)-Cd(1,1)*(T1-T2)));
EK12= Dt * ((1/C(1,2))*(Qb(1,2)-pC*BFB(1,2)*(T2-T25)+Cd(1,1)*(T1-T2)-Cd(1,2)*(T2-T3)));
EK13= Dt * ((1/C(1,3))*(Qb(1,3)-pC*BFB(1,3)*(T3-T25)+Cd(1,2)*(T2-T3)-Cd(1,3)*(T3-T4)));
EK14= Dt * ((1/C(1,4))*(Qb(1,4)-pC*((((BFB(1,4)+SKINV(1)*DL)/(1+SKINC(1)*ST))*2^(Err(4)/10))*(T4-T25))+Cd(1,3)*(T3-T4)-Qt(1)-E(1)));
B1=(pC*BFB(1,1)*(T1-T25))+(pC*BFB(1,2)*(T2-T25))+(BFB(1,3)*(T3-T25))+(pC*((((BFB(1,4)+SKINV(1)*DL)/(1+SKINC(1)*ST))*2^(Err(4)/10))*(T4-T25)));
% f(4,1) represents the Head Skin layer, with DL, ST, Err, Qt & E as variables %

%TRUNK SEGMENT% 
EK15= Dt * ((1/C(2,1))*(Qb(2,1)-pC*BFB(2,1)*(T5-T25)-Cd(2,1)*(T5-T6))+((0.0014*(34-ta)+ 0.017*(5.867-pa))*(sum(sum(Qb))))); %MISSING RES, heat loss by respiration%
EK16= Dt * ((1/C(2,2))*(Qb(2,2)-pC*BFB(2,2)*(T6-T25)+Cd(2,1)*(T5-T6)-Cd(2,2)*(T6-T7)));
EK17= Dt * ((1/C(2,3))*(Qb(2,3)-pC*BFB(2,3)*(T7-T25)+Cd(2,2)*(T6-T7)-Cd(2,3)*(T7-T8)));
EK18= Dt * ((1/C(2,4))*(Qb(2,4)-pC*((((BFB(2,4)+SKINV(2)*DL)/(1+SKINC(2)*ST))*2^(Err(8)/10))*(T8-T25))+Cd(2,3)*(T7-T8)-Qt(2)-E(2)));
B2=(pC*BFB(2,1)*(T5-T25))+(pC*BFB(2,2)*(T6-T25))+(pC*BFB(2,3)*(T7-T25))+(pC*((((BFB(2,4)+SKINV(2)*DL)/(1+SKINC(2)*ST))*2^(Err(8)/10))*(T8-T25)));

%ARMS SEGMENT%
EK19= Dt * ((1/C(3,1))*(Qb(3,1)-pC*BFB(3,1)*(T9-T25)-Cd(3,1)*(T9-T10)));
EK110= Dt * ((1/C(3,2))*(Qb(3,2)-pC*BFB(3,2)*(T10-T25)+Cd(3,1)*(T9-T10)-Cd(3,2)*(T10-T11)));
EK111= Dt * ((1/C(3,3))*(Qb(3,3)-pC*BFB(3,3)*(T11-T25)+Cd(3,2)*(T10-T11)-Cd(3,3)*(T11-T12)));
EK112= Dt * ((1/C(3,4))*(Qb(3,4)-pC*((((BFB(3,4)+SKINV(3)*DL)/(1+SKINC(3)*ST))*2^(Err(12)/10))*(T12-T25))+Cd(3,3)*(T11-T12)-Qt(3)-E(3)));
B3=(pC*BFB(3,1)*(T9-T25))+(pC*BFB(3,2)*(T10-T25))+(pC*BFB(3,3)*(T11-T25))+(pC*((((BFB(3,4)+SKINV(3)*DL)/(1+SKINC(3)*ST))*2^(Err(12)/10))*(T12-T25)));

%Hands%
EK113= Dt * ((1/C(4,1))*(Qb(4,1)-pC*BFB(4,1)*(T13-T25)-Cd(4,1)*(T13-T14)));
EK114= Dt * ((1/C(4,2))*(Qb(4,2)-pC*BFB(4,2)*(T14-T25)+Cd(4,1)*(T13-T14)-Cd(4,2)*(T14-T15)));
EK115= Dt * ((1/C(4,3))*(Qb(4,3)-pC*BFB(4,3)*(T15-T25)+Cd(4,2)*(T14-T15)-Cd(4,3)*(T15-T16)));
EK116= Dt * ((1/C(4,4))*(Qb(4,4)-pC*((((BFB(4,4)+SKINV(4)*DL)/(1+SKINC(4)*ST))*2^(Err(16)/10))*(T16-T25))+Cd(4,3)*(T15-T16)-Qt(4)-E(4)));
B4=(pC*BFB(4,1)*(T13-T25))+(pC*BFB(4,2)*(T14-T25))+(pC*BFB(4,3)*(T15-T25))+(pC*((((BFB(4,4)+SKINV(4)*DL)/(1+SKINC(4)*ST))*2^(Err(16)/10))*(T16-T25)));

%LEGS%
EK117= Dt * ((1/C(5,1))*(Qb(5,1)-pC*BFB(5,1)*(T17-T25)-Cd(5,1)*(T17-T18)));
EK118= Dt * ((1/C(5,2))*(Qb(5,2)-pC*BFB(5,2)*(T18-T25)+Cd(5,1)*(T17-T18)-Cd(5,2)*(T18-T19)));
EK119= Dt * ((1/C(5,3))*(Qb(5,3)-pC*BFB(5,3)*(T19-T25)+Cd(5,2)*(T18-T19)-Cd(5,3)*(T19-T20)));
EK120= Dt * ((1/C(5,4))*(Qb(5,4)-pC*((((BFB(5,4)+SKINV(5)*DL)/(1+SKINC(5)*ST))*2^(Err(20)/10))*(T20-T25))+Cd(5,3)*(T19-T20)-Qt(5)-E(5)));
B5=(pC*BFB(5,1)*(T17-T25))+(pC*BFB(5,2)*(T18-T25))+(pC*BFB(5,3)*(T19-T25))+(pC*((((BFB(5,4)+SKINV(5)*DL)/(1+SKINC(5)*ST))*2^(Err(20)/10))*(T20-T25)));

%FEET%
EK121= Dt * ((1/C(6,1))*(Qb(6,1)-pC*BFB(6,1)*(T21-T25)-Cd(6,1)*(T21-T22)));
EK122= Dt * ((1/C(6,2))*(Qb(6,2)-pC*BFB(6,2)*(T22-T25)+Cd(6,1)*(T21-T22)-Cd(6,2)*(T22-T23)));
EK123= Dt * ((1/C(6,3))*(Qb(6,3)-pC*BFB(6,3)*(T23-T25)+Cd(6,2)*(T22-T23)-Cd(6,3)*(T23-T24)));
EK124= Dt * ((1/C(6,4))*(Qb(6,4)-pC*((((BFB(6,4)+SKINV(6)*DL)/(1+SKINC(6)*ST))*2^(Err(24)/10))*(T24-T25))+Cd(6,3)*(T23-T24)-Qt(6)-E(6)));
B6=(pC*BFB(6,1)*(T21-T25))+(pC*BFB(6,2)*(T22-T25))+(pC*BFB(6,3)*(T23-T25))+(pC*((((BFB(6,4)+SKINV(6)*DL)/(1+SKINC(6)*ST))*2^(Err(24)/10))*(T24-T25)));

%Central blood ODE comes last%

EK125= Dt * (1/C(7,1))*(B1+B2+B3+B4+B5+B6);

T1=y1+EK11/2;
T2=y2+EK12/2;
T3=y3+EK13/2;
T4=y4+EK14/2;
T5=y5+EK15/2;
T6=y6+EK16/2;
T7=y7+EK17/2;
T8=y8+EK18/2;
T9=y9+EK19/2;
T10=y10+EK110/2;
T11=y11+EK111/2;
T12=y12+EK112/2;
T13=y13+EK113/2;
T14=y14+EK114/2;
T15=y15+EK115/2;
T16=y16+EK116/2;
T17=y17+EK117/2;
T18=y18+EK118/2;
T19=y19+EK119/2;
T20=y20+EK120/2;
T21=y21+EK121/2;
T22=y22+EK122/2;
T23=y23+EK123/2;
T24=y24+EK124/2;
T25=y25+EK125/2;

%HEAD SEGMENT%
EK21= Dt * ((1/C(1,1))*(Qb(1,1)-pC*BFB(1,1)*(T1-T25)-Cd(1,1)*(T1-T2)));
EK22= Dt * ((1/C(1,2))*(Qb(1,2)-pC*BFB(1,2)*(T2-T25)+Cd(1,1)*(T1-T2)-Cd(1,2)*(T2-T3)));
EK23= Dt * ((1/C(1,3))*(Qb(1,3)-pC*BFB(1,3)*(T3-T25)+Cd(1,2)*(T2-T3)-Cd(1,3)*(T3-T4)));
EK24= Dt * ((1/C(1,4))*(Qb(1,4)-pC*((((BFB(1,4)+SKINV(1)*DL)/(1+SKINC(1)*ST))*2^(Err(4)/10))*(T4-T25))+Cd(1,3)*(T3-T4)-Qt(1)-E(1)));
B7=(pC*BFB(1,1)*(T1-T25))+(pC*BFB(1,2)*(T2-T25))+(BFB(1,3)*(T3-T25))+(pC*((((BFB(1,4)+SKINV(1)*DL)/(1+SKINC(1)*ST))*2^(Err(4)/10))*(T4-T25)));

% f(4,1) represents the Head Skin layer, with DL, ST, Err, Qt & E as variables %

%TRUNK SEGMENT% 
EK25= Dt * ((1/C(2,1))*(Qb(2,1)-pC*BFB(2,1)*(T5-T25)-Cd(2,1)*(T5-T6))+((0.0014*(34-ta)+ 0.017*(5.867-pa))*(sum(sum(Qb))))); %MISSING RES, heat loss by respiration%
EK26= Dt * ((1/C(2,2))*(Qb(2,2)-pC*BFB(2,2)*(T6-T25)+Cd(2,1)*(T5-T6)-Cd(2,2)*(T6-T7)));
EK27= Dt * ((1/C(2,3))*(Qb(2,3)-pC*BFB(2,3)*(T7-T25)+Cd(2,2)*(T6-T7)-Cd(2,3)*(T7-T8)));
EK28= Dt * ((1/C(2,4))*(Qb(2,4)-pC*((((BFB(2,4)+SKINV(2)*DL)/(1+SKINC(2)*ST))*2^(Err(8)/10))*(T8-T25))+Cd(2,3)*(T7-T8)-Qt(2)-E(2)));
B8=(pC*BFB(2,1)*(T5-T25))+(pC*BFB(2,2)*(T6-T25))+(pC*BFB(2,3)*(T7-T25))+(pC*((((BFB(2,4)+SKINV(2)*DL)/(1+SKINC(2)*ST))*2^(Err(8)/10))*(T8-T25)));

%ARMS SEGMENT%
EK29= Dt * ((1/C(3,1))*(Qb(3,1)-pC*BFB(3,1)*(T9-T25)-Cd(3,1)*(T9-T10)));
EK210= Dt * ((1/C(3,2))*(Qb(3,2)-pC*BFB(3,2)*(T10-T25)+Cd(3,1)*(T9-T10)-Cd(3,2)*(T10-T11)));
EK211= Dt * ((1/C(3,3))*(Qb(3,3)-pC*BFB(3,3)*(T11-T25)+Cd(3,2)*(T10-T11)-Cd(3,3)*(T11-T12)));
EK212= Dt * ((1/C(3,4))*(Qb(3,4)-pC*((((BFB(3,4)+SKINV(3)*DL)/(1+SKINC(3)*ST))*2^(Err(12)/10))*(T12-T25))+Cd(3,3)*(T11-T12)-Qt(3)-E(3)));
B9=(pC*BFB(3,1)*(T9-T25))+(pC*BFB(3,2)*(T10-T25))+(pC*BFB(3,3)*(T11-T25))+(pC*((((BFB(3,4)+SKINV(3)*DL)/(1+SKINC(3)*ST))*2^(Err(12)/10))*(T12-T25)));

%Hands%
EK213= Dt * ((1/C(4,1))*(Qb(4,1)-pC*BFB(4,1)*(T13-T25)-Cd(4,1)*(T13-T14)));
EK214= Dt * ((1/C(4,2))*(Qb(4,2)-pC*BFB(4,2)*(T14-T25)+Cd(4,1)*(T13-T14)-Cd(4,2)*(T14-T15)));
EK215= Dt * ((1/C(4,3))*(Qb(4,3)-pC*BFB(4,3)*(T15-T25)+Cd(4,2)*(T14-T15)-Cd(4,3)*(T15-T16)));
EK216= Dt * ((1/C(4,4))*(Qb(4,4)-pC*((((BFB(4,4)+SKINV(4)*DL)/(1+SKINC(4)*ST))*2^(Err(16)/10))*(T16-T25))+Cd(4,3)*(T15-T16)-Qt(4)-E(4)));
B10=(pC*BFB(4,1)*(T13-T25))+(pC*BFB(4,2)*(T14-T25))+(pC*BFB(4,3)*(T15-T25))+(pC*((((BFB(4,4)+SKINV(4)*DL)/(1+SKINC(4)*ST))*2^(Err(16)/10))*(T16-T25)));

%LEGS%
EK217= Dt * ((1/C(5,1))*(Qb(5,1)-pC*BFB(5,1)*(T17-T25)-Cd(5,1)*(T17-T18)));
EK218= Dt * ((1/C(5,2))*(Qb(5,2)-pC*BFB(5,2)*(T18-T25)+Cd(5,1)*(T17-T18)-Cd(5,2)*(T18-T19)));
EK219= Dt * ((1/C(5,3))*(Qb(5,3)-pC*BFB(5,3)*(T19-T25)+Cd(5,2)*(T18-T19)-Cd(5,3)*(T19-T20)));
EK220= Dt * ((1/C(5,4))*(Qb(5,4)-pC*((((BFB(5,4)+SKINV(5)*DL)/(1+SKINC(5)*ST))*2^(Err(20)/10))*(T20-T25))+Cd(5,3)*(T19-T20)-Qt(5)-E(5)));
B11=(pC*BFB(5,1)*(T17-T25))+(pC*BFB(5,2)*(T18-T25))+(pC*BFB(5,3)*(T19-T25))+(pC*((((BFB(5,4)+SKINV(5)*DL)/(1+SKINC(5)*ST))*2^(Err(20)/10))*(T20-T25)));

%FEET%
EK221= Dt * ((1/C(6,1))*(Qb(6,1)-pC*BFB(6,1)*(T21-T25)-Cd(6,1)*(T21-T22)));
EK222= Dt * ((1/C(6,2))*(Qb(6,2)-pC*BFB(6,2)*(T22-T25)+Cd(6,1)*(T21-T22)-Cd(6,2)*(T22-T23)));
EK223= Dt * ((1/C(6,3))*(Qb(6,3)-pC*BFB(6,3)*(T23-T25)+Cd(6,2)*(T22-T23)-Cd(6,3)*(T23-T24)));
EK224= Dt * ((1/C(6,4))*(Qb(6,4)-pC*((((BFB(6,4)+SKINV(6)*DL)/(1+SKINC(6)*ST))*2^(Err(24)/10))*(T24-T25))+Cd(6,3)*(T23-T24)-Qt(6)-E(6)));
B12=(pC*BFB(6,1)*(T21-T25))+(pC*BFB(6,2)*(T22-T25))+(pC*BFB(6,3)*(T23-T25))+(pC*((((BFB(6,4)+SKINV(6)*DL)/(1+SKINC(6)*ST))*2^(Err(24)/10))*(T24-T25)));

%Central blood ODE comes last%

EK225= Dt * (1/C(7,1))*(B7+B8+B9+B10+B11+B12);

T1=y1+EK21/2;
T2=y2+EK22/2;
T3=y3+EK23/2;
T4=y4+EK24/2;
T5=y5+EK25/2;
T6=y6+EK26/2;
T7=y7+EK27/2;
T8=y8+EK28/2;
T9=y9+EK29/2;
T10=y10+EK210/2;
T11=y11+EK211/2;
T12=y12+EK212/2;
T13=y13+EK213/2;
T14=y14+EK214/2;
T15=y15+EK215/2;
T16=y16+EK216/2;
T17=y17+EK217/2;
T18=y18+EK218/2;
T19=y19+EK219/2;
T20=y20+EK220/2;
T21=y21+EK221/2;
T22=y22+EK222/2;
T23=y23+EK223/2;
T24=y24+EK224/2;
T25=y25+EK225/2;

%HEAD SEGMENT%
EK31= Dt * ((1/C(1,1))*(Qb(1,1)-pC*BFB(1,1)*(T1-T25)-Cd(1,1)*(T1-T2)));
EK32= Dt * ((1/C(1,2))*(Qb(1,2)-pC*BFB(1,2)*(T2-T25)+Cd(1,1)*(T1-T2)-Cd(1,2)*(T2-T3)));
EK33= Dt * ((1/C(1,3))*(Qb(1,3)-pC*BFB(1,3)*(T3-T25)+Cd(1,2)*(T2-T3)-Cd(1,3)*(T3-T4)));
EK34= Dt * ((1/C(1,4))*(Qb(1,4)-pC*((((BFB(1,4)+SKINV(1)*DL)/(1+SKINC(1)*ST))*2^(Err(4)/10))*(T4-T25))+Cd(1,3)*(T3-T4)-Qt(1)-E(1)));
B13=(pC*BFB(1,1)*(T1-T25))+(pC*BFB(1,2)*(T2-T25))+(BFB(1,3)*(T3-T25))+(pC*((((BFB(1,4)+SKINV(1)*DL)/(1+SKINC(1)*ST))*2^(Err(4)/10))*(T4-T25)));

% f(4,1) represents the Head Skin layer, with DL, ST, Err, Qt & E as variables %

%TRUNK SEGMENT% 
EK35= Dt * ((1/C(2,1))*(Qb(2,1)-pC*BFB(2,1)*(T5-T25)-Cd(2,1)*(T5-T6))+((0.0014*(34-ta)+ 0.017*(5.867-pa))*(sum(sum(Qb))))); %MISSING RES, heat loss by respiration%
EK36= Dt * ((1/C(2,2))*(Qb(2,2)-pC*BFB(2,2)*(T6-T25)+Cd(2,1)*(T5-T6)-Cd(2,2)*(T6-T7)));
EK37= Dt * ((1/C(2,3))*(Qb(2,3)-pC*BFB(2,3)*(T7-T25)+Cd(2,2)*(T6-T7)-Cd(2,3)*(T7-T8)));
EK38= Dt * ((1/C(2,4))*(Qb(2,4)-pC*((((BFB(2,4)+SKINV(2)*DL)/(1+SKINC(2)*ST))*2^(Err(8)/10))*(T8-T25))+Cd(2,3)*(T7-T8)-Qt(2)-E(2)));
B14=(pC*BFB(2,1)*(T5-T25))+(pC*BFB(2,2)*(T6-T25))+(pC*BFB(2,3)*(T7-T25))+(pC*((((BFB(2,4)+SKINV(2)*DL)/(1+SKINC(2)*ST))*2^(Err(8)/10))*(T8-T25)));

%ARMS SEGMENT%
EK39= Dt * ((1/C(3,1))*(Qb(3,1)-pC*BFB(3,1)*(T9-T25)-Cd(3,1)*(T9-T10)));
EK310= Dt * ((1/C(3,2))*(Qb(3,2)-pC*BFB(3,2)*(T10-T25)+Cd(3,1)*(T9-T10)-Cd(3,2)*(T10-T11)));
EK311= Dt * ((1/C(3,3))*(Qb(3,3)-pC*BFB(3,3)*(T11-T25)+Cd(3,2)*(T10-T11)-Cd(3,3)*(T11-T12)));
EK312= Dt * ((1/C(3,4))*(Qb(3,4)-pC*((((BFB(3,4)+SKINV(3)*DL)/(1+SKINC(3)*ST))*2^(Err(12)/10))*(T12-T25))+Cd(3,3)*(T11-T12)-Qt(3)-E(3)));
B15=(pC*BFB(3,1)*(T9-T25))+(pC*BFB(3,2)*(T10-T25))+(pC*BFB(3,3)*(T11-T25))+(pC*((((BFB(3,4)+SKINV(3)*DL)/(1+SKINC(3)*ST))*2^(Err(12)/10))*(T12-T25)));

%Hands%
EK313= Dt * ((1/C(4,1))*(Qb(4,1)-pC*BFB(4,1)*(T13-T25)-Cd(4,1)*(T13-T14)));
EK314= Dt * ((1/C(4,2))*(Qb(4,2)-pC*BFB(4,2)*(T14-T25)+Cd(4,1)*(T13-T14)-Cd(4,2)*(T14-T15)));
EK315= Dt * ((1/C(4,3))*(Qb(4,3)-pC*BFB(4,3)*(T15-T25)+Cd(4,2)*(T14-T15)-Cd(4,3)*(T15-T16)));
EK316= Dt * ((1/C(4,4))*(Qb(4,4)-pC*((((BFB(4,4)+SKINV(4)*DL)/(1+SKINC(4)*ST))*2^(Err(16)/10))*(T16-T25))+Cd(4,3)*(T15-T16)-Qt(4)-E(4)));
B16=(pC*BFB(4,1)*(T13-T25))+(pC*BFB(4,2)*(T14-T25))+(pC*BFB(4,3)*(T15-T25))+(pC*((((BFB(4,4)+SKINV(4)*DL)/(1+SKINC(4)*ST))*2^(Err(16)/10))*(T16-T25)));

%LEGS%
EK317= Dt * ((1/C(5,1))*(Qb(5,1)-pC*BFB(5,1)*(T17-T25)-Cd(5,1)*(T17-T18)));
EK318= Dt * ((1/C(5,2))*(Qb(5,2)-pC*BFB(5,2)*(T18-T25)+Cd(5,1)*(T17-T18)-Cd(5,2)*(T18-T19)));
EK319= Dt * ((1/C(5,3))*(Qb(5,3)-pC*BFB(5,3)*(T19-T25)+Cd(5,2)*(T18-T19)-Cd(5,3)*(T19-T20)));
EK320= Dt * ((1/C(5,4))*(Qb(5,4)-pC*((((BFB(5,4)+SKINV(5)*DL)/(1+SKINC(5)*ST))*2^(Err(20)/10))*(T20-T25))+Cd(5,3)*(T19-T20)-Qt(5)-E(5)));
B17=(pC*BFB(5,1)*(T17-T25))+(pC*BFB(5,2)*(T18-T25))+(pC*BFB(5,3)*(T19-T25))+(pC*((((BFB(5,4)+SKINV(5)*DL)/(1+SKINC(5)*ST))*2^(Err(20)/10))*(T20-T25)));

%FEET%
EK321= Dt * ((1/C(6,1))*(Qb(6,1)-pC*BFB(6,1)*(T21-T25)-Cd(6,1)*(T21-T22)));
EK322= Dt * ((1/C(6,2))*(Qb(6,2)-pC*BFB(6,2)*(T22-T25)+Cd(6,1)*(T21-T22)-Cd(6,2)*(T22-T23)));
EK323= Dt * ((1/C(6,3))*(Qb(6,3)-pC*BFB(6,3)*(T23-T25)+Cd(6,2)*(T22-T23)-Cd(6,3)*(T23-T24)));
EK324= Dt * ((1/C(6,4))*(Qb(6,4)-pC*((((BFB(6,4)+SKINV(6)*DL)/(1+SKINC(6)*ST))*2^(Err(24)/10))*(T24-T25))+Cd(6,3)*(T23-T24)-Qt(6)-E(6)));
B18=(pC*BFB(6,1)*(T21-T25))+(pC*BFB(6,2)*(T22-T25))+(pC*BFB(6,3)*(T23-T25))+(pC*((((BFB(6,4)+SKINV(6)*DL)/(1+SKINC(6)*ST))*2^(Err(24)/10))*(T24-T25)));

%Central blood ODE comes last%

EK325= Dt * (1/C(7,1))*(B13+B14+B15+B16+B17+B18);

T1=y1+EK31;
T2=y2+EK32;
T3=y3+EK33;
T4=y4+EK34;
T5=y5+EK35;
T6=y6+EK36;
T7=y7+EK37;
T8=y8+EK38;
T9=y9+EK39;
T10=y10+EK310;
T11=y11+EK311;
T12=y12+EK312;
T13=y13+EK313;
T14=y14+EK314;
T15=y15+EK315;
T16=y16+EK316;
T17=y17+EK317;
T18=y18+EK318;
T19=y19+EK319;
T20=y20+EK320;
T21=y21+EK321;
T22=y22+EK322;
T23=y23+EK323;
T24=y24+EK324;
T25=y25+EK325;

%HEAD SEGMENT%
EK41= Dt * ((1/C(1,1))*(Qb(1,1)-pC*BFB(1,1)*(T1-T25)-Cd(1,1)*(T1-T2)));
EK42= Dt * ((1/C(1,2))*(Qb(1,2)-pC*BFB(1,2)*(T2-T25)+Cd(1,1)*(T1-T2)-Cd(1,2)*(T2-T3)));
EK43= Dt * ((1/C(1,3))*(Qb(1,3)-pC*BFB(1,3)*(T3-T25)+Cd(1,2)*(T2-T3)-Cd(1,3)*(T3-T4)));
EK44= Dt * ((1/C(1,4))*(Qb(1,4)-pC*((((BFB(1,4)+SKINV(1)*DL)/(1+SKINC(1)*ST))*2^(Err(4)/10))*(T4-T25))+Cd(1,3)*(T3-T4)-Qt(1)-E(1)));
%sum of B
B19=(pC*BFB(1,1)*(T1-T25))+(pC*BFB(1,2)*(T2-T25))+(BFB(1,3)*(T3-T25))+(pC*((((BFB(1,4)+SKINV(1)*DL)/(1+SKINC(1)*ST))*2^(Err(4)/10))*(T4-T25)));

% f(4,1) represents the Head Skin layer, with DL, ST, Err, Qt & E as variables %

%TRUNK SEGMENT% 
EK45= Dt * ((1/C(2,1))*(Qb(2,1)-pC*BFB(2,1)*(T5-T25)-Cd(2,1)*(T5-T6))+((0.0014*(34-ta)+ 0.017*(5.867-pa))*(sum(sum(Qb))))); %MISSING RES, heat loss by respiration%
EK46= Dt * ((1/C(2,2))*(Qb(2,2)-pC*BFB(2,2)*(T6-T25)+Cd(2,1)*(T5-T6)-Cd(2,2)*(T6-T7)));
EK47= Dt * ((1/C(2,3))*(Qb(2,3)-pC*BFB(2,3)*(T7-T25)+Cd(2,2)*(T6-T7)-Cd(2,3)*(T7-T8)));
EK48= Dt * ((1/C(2,4))*(Qb(2,4)-pC*((((BFB(2,4)+SKINV(2)*DL)/(1+SKINC(2)*ST))*2^(Err(8)/10))*(T8-T25))+Cd(2,3)*(T7-T8)-Qt(2)-E(2)));
B20=(pC*BFB(2,1)*(T5-T25))+(pC*BFB(2,2)*(T6-T25))+(pC*BFB(2,3)*(T7-T25))+(pC*((((BFB(2,4)+SKINV(2)*DL)/(1+SKINC(2)*ST))*2^(Err(8)/10))*(T8-T25)));

%ARMS SEGMENT%
EK49= Dt * ((1/C(3,1))*(Qb(3,1)-pC*BFB(3,1)*(T9-T25)-Cd(3,1)*(T9-T10)));
EK410= Dt * ((1/C(3,2))*(Qb(3,2)-pC*BFB(3,2)*(T10-T25)+Cd(3,1)*(T9-T10)-Cd(3,2)*(T10-T11)));
EK411= Dt * ((1/C(3,3))*(Qb(3,3)-pC*BFB(3,3)*(T11-T25)+Cd(3,2)*(T10-T11)-Cd(3,3)*(T11-T12)));
EK412= Dt * ((1/C(3,4))*(Qb(3,4)-pC*((((BFB(3,4)+SKINV(3)*DL)/(1+SKINC(3)*ST))*2^(Err(12)/10))*(T12-T25))+Cd(3,3)*(T11-T12)-Qt(3)-E(3)));
B21=(pC*BFB(3,1)*(T9-T25))+(pC*BFB(3,2)*(T10-T25))+(pC*BFB(3,3)*(T11-T25))+(pC*((((BFB(3,4)+SKINV(3)*DL)/(1+SKINC(3)*ST))*2^(Err(12)/10))*(T12-T25)));

%Hands%
EK413= Dt * ((1/C(4,1))*(Qb(4,1)-pC*BFB(4,1)*(T13-T25)-Cd(4,1)*(T13-T14)));
EK414= Dt * ((1/C(4,2))*(Qb(4,2)-pC*BFB(4,2)*(T14-T25)+Cd(4,1)*(T13-T14)-Cd(4,2)*(T14-T15)));
EK415= Dt * ((1/C(4,3))*(Qb(4,3)-pC*BFB(4,3)*(T15-T25)+Cd(4,2)*(T14-T15)-Cd(4,3)*(T15-T16)));
EK416= Dt * ((1/C(4,4))*(Qb(4,4)-pC*((((BFB(4,4)+SKINV(4)*DL)/(1+SKINC(4)*ST))*2^(Err(16)/10))*(T16-T25))+Cd(4,3)*(T15-T16)-Qt(4)-E(4)));
B22=(pC*BFB(4,1)*(T13-T25))+(pC*BFB(4,2)*(T14-T25))+(pC*BFB(4,3)*(T15-T25))+(pC*((((BFB(4,4)+SKINV(4)*DL)/(1+SKINC(4)*ST))*2^(Err(16)/10))*(T16-T25)));

%LEGS%
EK417= Dt * ((1/C(5,1))*(Qb(5,1)-pC*BFB(5,1)*(T17-T25)-Cd(5,1)*(T17-T18)));
EK418= Dt * ((1/C(5,2))*(Qb(5,2)-pC*BFB(5,2)*(T18-T25)+Cd(5,1)*(T17-T18)-Cd(5,2)*(T18-T19)));
EK419= Dt * ((1/C(5,3))*(Qb(5,3)-pC*BFB(5,3)*(T19-T25)+Cd(5,2)*(T18-T19)-Cd(5,3)*(T19-T20)));
EK420= Dt * ((1/C(5,4))*(Qb(5,4)-pC*((((BFB(5,4)+SKINV(5)*DL)/(1+SKINC(5)*ST))*2^(Err(20)/10))*(T20-T25))+Cd(5,3)*(T19-T20)-Qt(5)-E(5)));
B23=(pC*BFB(5,1)*(T17-T25))+(pC*BFB(5,2)*(T18-T25))+(pC*BFB(5,3)*(T19-T25))+(pC*((((BFB(5,4)+SKINV(5)*DL)/(1+SKINC(5)*ST))*2^(Err(20)/10))*(T20-T25)));

%FEET%
EK421= Dt * ((1/C(6,1))*(Qb(6,1)-pC*BFB(6,1)*(T21-T25)-Cd(6,1)*(T21-T22)));
EK422= Dt * ((1/C(6,2))*(Qb(6,2)-pC*BFB(6,2)*(T22-T25)+Cd(6,1)*(T21-T22)-Cd(6,2)*(T22-T23)));
EK423= Dt * ((1/C(6,3))*(Qb(6,3)-pC*BFB(6,3)*(T23-T25)+Cd(6,2)*(T22-T23)-Cd(6,3)*(T23-T24)));
EK424= Dt * ((1/C(6,4))*(Qb(6,4)-pC*((((BFB(6,4)+SKINV(6)*DL)/(1+SKINC(6)*ST))*2^(Err(24)/10))*(T24-T25))+Cd(6,3)*(T23-T24)-Qt(6)-E(6)));
B24=(pC*BFB(6,1)*(T21-T25))+(pC*BFB(6,2)*(T22-T25))+(pC*BFB(6,3)*(T23-T25))+(pC*((((BFB(6,4)+SKINV(6)*DL)/(1+SKINC(6)*ST))*2^(Err(24)/10))*(T24-T25)));

%Central blood ODE comes last%

EK425= Dt * (1/C(7,1))*(B19+B20+B21+B22+B23+B24);

y1=y1+(EK11+2*EK21+2*EK31+EK41)/6;
y2=y2+(EK12+2*EK22+2*EK32+EK42)/6;
y3=y3+(EK13+2*EK23+2*EK33+EK43)/6;
y4=y4+(EK14+2*EK24+2*EK34+EK44)/6;
y5=y5+(EK15+2*EK25+2*EK35+EK45)/6;
y6=y6+(EK16+2*EK26+2*EK36+EK46)/6;
y7=y7+(EK17+2*EK27+2*EK37+EK47)/6;
y8=y8+(EK18+2*EK28+2*EK38+EK48)/6;
y9=y9+(EK19+2*EK29+2*EK39+EK49)/6;
y10=y10+(EK110+2*EK210+2*EK310+EK410)/6;
y11=y11+(EK111+2*EK211+2*EK311+EK411)/6;
y12=y12+(EK112+2*EK212+2*EK312+EK412)/6;
y13=y13+(EK113+2*EK213+2*EK313+EK413)/6;
y14=y14+(EK114+2*EK214+2*EK314+EK414)/6;
y15=y15+(EK115+2*EK215+2*EK315+EK415)/6;
y16=y16+(EK116+2*EK216+2*EK316+EK416)/6;
y17=y17+(EK117+2*EK217+2*EK317+EK417)/6;
y18=y18+(EK118+2*EK218+2*EK318+EK418)/6;
y19=y19+(EK119+2*EK219+2*EK319+EK419)/6;
y20=y20+(EK120+2*EK220+2*EK320+EK420)/6;
y21=y21+(EK121+2*EK221+2*EK321+EK421)/6;
y22=y22+(EK122+2*EK222+2*EK322+EK422)/6;
y23=y23+(EK123+2*EK223+2*EK323+EK423)/6;
y24=y24+(EK124+2*EK224+2*EK324+EK424)/6;
y25=y25+(EK125+2*EK225+2*EK325+EK425)/6;

%Final Values
T1=y1;
T2=y2;
T3=y3;
T4=y4;
T5=y5;
T6=y6;
T7=y7;
T8=y8;
T9=y9;
T10=y10;
T11=y11;
T12=y12;
T13=y13;
T14=y14;
T15=y15;
T16=y16;
T17=y17;
T18=y18;
T19=y19;
T20=y20;
T21=y21;
T22=y22;
T23=y23;
T24=y24;
T25=y25;

%Assigning to Matrix
MT1(i+1)=T1;
MT2(i+1)=T2;
MT3(i+1)=T3;
MT4(i+1)=T4;
MT5(i+1)=T5;
MT6(i+1)=T6;
MT7(i+1)=T7;
MT8(i+1)=T8;
MT9(i+1)=T9;
MT10(i+1)=T10;
MT11(i+1)=T11;
MT12(i+1)=T12;
MT13(i+1)=T13;
MT14(i+1)=T14;
MT15(i+1)=T15;
MT16(i+1)=T16;
MT17(i+1)=T17;
MT18(i+1)=T18;
MT19(i+1)=T19;
MT20(i+1)=T20;
MT21(i+1)=T21;
MT22(i+1)=T22;
MT23(i+1)=T23;
MT24(i+1)=T24;
MT25(i+1)=T25;

t(i+1)=a+i*Dt;

end
plot(t,MT1)

end