Finding difficult in solving these coupled seven ODE's.. Initially tried to solve by runge kutta method... I want to solve by through ODE45 solver because solving by rk method requires more lines of execute,

clear;

%Runge kutta for four coupled system of equations%

%depsi_dt=A/(1-omega)^(b*sigma(1-hc)*sinh[(B*sigma*(1-H))]

%dh_dt=(h*sigma)*(1-(H/Hstar)*(A/(1-omega)*sinh[(B*sigma*(1-H))]

%domega_dt=D*A/(1-omega)*sinh[(B*sigma*(1-H))]

%Define function handles

%Define function handles

ftf=@(t,tf,tg1,tg2,tag,tab,tpot,tload) (1.4915-1.077E-10*(tf^4-8.03E+8)-0.0408*(tf-305)+9.7E-11*(tf^4-tg1^4));

ftg1=@(t,tf,tg1,tg2,tag,tab,tpot,tload) (6.54E-3-3.2E-12*(tg1^4-8.03E+8)-1.5E-3*(tg1-305)-4.5E-14*(tf^4-tg1^4)+((1.26E-3*(tg2-tg1)^1.27)/((tg1+tg2)^0.31))+5.88E-12*(tg2^4-tg1^4));

ftg2=@(t,tf,tg1,tg2,tag,tab,tpot,tload) (7.104E-3+((0.295*((tag-tg2)^1.304)/((tag+tg2)^0.304))+2.087E-8*(tab^4-tg2^4)+0.261E-8*(tpot^4-tg2^4)-((1.26E-3*(tg2-tg1)^1.27))/((tg1+tg2)^0.31))-5.88E-12*(tg2^4-tg1^4));

ftag=@(t,tf,tg1,tg2,tag,tab,tpot,tload) (0.06*(tab-tag)^1.25-((1.3169*(tag-tg2)^1.304)/((tag+tg2)^0.304)))-6.37E-3*(tag-tpot)-((0.32*(tag-tpot)^1.166)/((tag+tpot)^0.166))-((1.12E-4*(tab-tag)*(tag-305))/(0.8703*(tab-tag)+0.606));

ftab=@(t,tf,tg1,tg2,tag,tab,tpot,tload) (-2.12E-11*(tab^4-tg2^4)-6.16E-4*(tab-tag)+8.103E-3-4.622E-4*(tab-305));

ftpot=@(t,tf,tg1,tg2,tag,tab,tpot,tload) (-3.536E-10*(tpot^4-tg2^4)+38.26*(tab-tpot)+8.635E-4*(tag-tpot)+((0.0442*(tag-tpot)^1.166)/((tag+tg2)^0.166))-((17.168*(tpot-tload)^1.25)/((tpot+tload)^0.25))+0.1646);

ftload=@(t,tf,tg1,tg2,tag,tab,tpot,tload) ((1.1835*(tpot-tload)^1.25)/((tpot+tload)^0.25));

%Initial conditions

t(1)=0;

tf(1)=305;

tg1(1)=304;

tg2(1)=305;

tag(1)=306;

tab(1)=306;

tpot(1)=304;

tload(1)=304;

%step size

h=10;

s=h/2;

tfinal=40;

N=ceil(tfinal/h);

%update loop

for i=1:N-1

%update time

t(i+1)=t(i)+h;

% Runge Kutta Method

K1tf=h*ftf(t(i) ,tf(i) ,tg1(i) ,tg2(i) ,tag(i) ,tab(i) ,tpot(i) ,tload(i));

K1tg1=h*ftg1(t(i) ,tf(i) ,tg1(i) ,tg2(i) ,tag(i) ,tab(i) ,tpot(i) ,tload(i));

K1tg2=h*ftg2(t(i) ,tf(i) ,tg1(i) ,tg2(i) ,tag(i) ,tab(i) ,tpot(i) ,tload(i));

K1tag=h*ftag(t(i) ,tf(i) ,tg1(i) ,tg2(i) ,tag(i) ,tab(i) ,tpot(i) ,tload(i));

K1tab=h*ftab(t(i) ,tf(i) ,tg1(i) ,tg2(i) ,tag(i) ,tab(i) ,tpot(i) ,tload(i));

K1tpot=h*ftpot(t(i) ,tf(i) ,tg1(i) ,tg2(i) ,tag(i) ,tab(i) ,tpot(i) ,tload(i));

K1tload=h*ftload(t(i) ,tf(i) ,tg1(i) ,tg2(i) ,tag(i) ,tab(i) ,tpot(i) ,tload(i));

K2tf=h*ftf(t(i)+s ,tf(i)+1/2*K1tf ,tg1(i)+1/2*K1tg1 ,tg2(i)+1/2*K1tg2 ,tag(i)+1/2*K1tag ,tab(i)+1/2*K1tab ,tpot(i)+1/2*K1tpot ,tload(i)+1/2*K1tload);

K2tg1=h*ftg1(t(i)+s ,tf(i)+1/2*K1tf ,tg1(i)+1/2*K1tg1 ,tg2(i)+1/2*K1tg2 ,tag(i)+1/2*K1tag ,tab(i)+1/2*K1tab ,tpot(i)+1/2*K1tpot ,tload(i)+1/2*K1tload);

K2tg2=h*ftg2(t(i)+s ,tf(i)+1/2*K1tf ,tg1(i)+1/2*K1tg1 ,tg2(i)+1/2*K1tg2 ,tag(i)+1/2*K1tag ,tab(i)+1/2*K1tab ,tpot(i)+1/2*K1tpot ,tload(i)+1/2*K1tload);

K2tag=h*ftag(t(i)+s ,tf(i)+1/2*K1tf ,tg1(i)+1/2*K1tg1 ,tg2(i)+1/2*K1tg2 ,tag(i)+1/2*K1tag ,tab(i)+1/2*K1tab ,tpot(i)+1/2*K1tpot ,tload(i)+1/2*K1tload);

K2tab=h*ftab(t(i)+s ,tf(i)+1/2*K1tf ,tg1(i)+1/2*K1tg1 ,tg2(i)+1/2*K1tg2 ,tag(i)+1/2*K1tag ,tab(i)+1/2*K1tab ,tpot(i)+1/2*K1tpot ,tload(i)+1/2*K1tload);

K2tpot=h*ftpot(t(i)+s ,tf(i)+1/2*K1tf ,tg1(i)+1/2*K1tg1 ,tg2(i)+1/2*K1tg2 ,tag(i)+1/2*K1tag ,tab(i)+1/2*K1tab ,tpot(i)+1/2*K1tpot ,tload(i)+1/2*K1tload);

K2tload=h*ftload(t(i)+s ,tf(i)+1/2*K1tf ,tg1(i)+1/2*K1tg1 ,tg2(i)+1/2*K1tg2 ,tag(i)+1/2*K1tag ,tab(i)+1/2*K1tab ,tpot(i)+1/2*K1tpot ,tload(i)+1/2*K1tload);

K3tf=h*ftf(t(i)+s ,tf(i)+1/2*K2tf ,tg1(i)+1/2*K2tg1 ,tg2(i)+1/2*K2tg2 ,tag(i)+1/2*K2tag ,tab(i)+1/2*K2tab ,tpot(i)+1/2*K2tpot ,tload(i)+1/2*K2tload);

K3tg1=h*ftg1(t(i)+s ,tf(i)+1/2*K2tf ,tg1(i)+1/2*K2tg1 ,tg2(i)+1/2*K2tg2 ,tag(i)+1/2*K2tag ,tab(i)+1/2*K2tab ,tpot(i)+1/2*K2tpot ,tload(i)+1/2*K2tload);

K3tg2=h*ftg2(t(i)+s ,tf(i)+1/2*K2tf ,tg1(i)+1/2*K2tg1 ,tg2(i)+1/2*K2tg2 ,tag(i)+1/2*K2tag ,tab(i)+1/2*K2tab ,tpot(i)+1/2*K2tpot ,tload(i)+1/2*K2tload);

K3tag=h*ftag(t(i)+s ,tf(i)+1/2*K2tf ,tg1(i)+1/2*K2tg1 ,tg2(i)+1/2*K2tg2 ,tag(i)+1/2*K2tag ,tab(i)+1/2*K2tab ,tpot(i)+1/2*K2tpot ,tload(i)+1/2*K2tload);

K3tab=h*ftab(t(i)+s ,tf(i)+1/2*K2tf ,tg1(i)+1/2*K2tg1 ,tg2(i)+1/2*K2tg2 ,tag(i)+1/2*K2tag ,tab(i)+1/2*K2tab ,tpot(i)+1/2*K2tpot ,tload(i)+1/2*K2tload);

K3tpot=h*ftpot(t(i)+s ,tf(i)+1/2*K2tf ,tg1(i)+1/2*K2tg1 ,tg2(i)+1/2*K2tg2 ,tag(i)+1/2*K2tag ,tab(i)+1/2*K2tab ,tpot(i)+1/2*K2tpot ,tload(i)+1/2*K2tload);

K3tload=h*ftload(t(i)+s ,tf(i)+1/2*K2tf ,tg1(i)+1/2*K2tg1 ,tg2(i)+1/2*K2tg2 ,tag(i)+1/2*K2tag ,tab(i)+1/2*K2tab ,tpot(i)+1/2*K2tpot ,tload(i)+1/2*K2tload);

K4tf=h*ftf(t(i)+h ,tf(i)+h*K3tf ,tg1(i)+h*K3tg1 ,tg2(i)+h*K3tg2 ,tag(i)+h*K3tag ,tab(i)+h*K3tab ,tpot(i)+h*K3tpot ,tload(i)+h*K3tload);

K4tg1=h*ftg1(t(i)+h ,tf(i)+h*K3tf ,tg1(i)+h*K3tg1 ,tg2(i)+h*K3tg2 ,tag(i)+h*K3tag ,tab(i)+h*K3tab ,tpot(i)+h*K3tpot ,tload(i)+h*K3tload);

K4tg2=h*ftg2(t(i)+h ,tf(i)+h*K3tf ,tg1(i)+h*K3tg1 ,tg2(i)+h*K3tg2 ,tag(i)+h*K3tag ,tab(i)+h*K3tab ,tpot(i)+h*K3tpot ,tload(i)+h*K3tload);

K4tag=h*ftag(t(i)+h ,tf(i)+h*K3tf ,tg1(i)+h*K3tg1 ,tg2(i)+h*K3tg2 ,tag(i)+h*K3tag ,tab(i)+h*K3tab ,tpot(i)+h*K3tpot ,tload(i)+h*K3tload);

K4tab=h*ftab(t(i)+h ,tf(i)+h*K3tf ,tg1(i)+h*K3tg1 ,tg2(i)+h*K3tg2 ,tag(i)+h*K3tag ,tab(i)+h*K3tab ,tpot(i)+h*K3tpot ,tload(i)+h*K3tload);

K4tpot=h*ftpot(t(i)+h ,tf(i)+h*K3tf ,tg1(i)+h*K3tg1 ,tg2(i)+h*K3tg2 ,tag(i)+h*K3tag ,tab(i)+h*K3tab ,tpot(i)+h*K3tpot ,tload(i)+h*K3tload);

K4tload=h*ftload(t(i)+h ,tf(i)+h*K3tf ,tg1(i)+h*K3tg1 ,tg2(i)+h*K3tg2 ,tag(i)+h*K3tag ,tab(i)+h*K3tab ,tpot(i)+h*K3tpot ,tload(i)+h*K3tload);

%% Repeating iteration step by step

tf(i+1)=tf(i)+1/6*(K1tf+2*K2tf+2*K3tf+K4tf);

tg1(i+1)=tg1(i)+1/6*(K1tg1+2*K2tg1+2*K3tg1+K4tg1);

tg2(i+1)=tg2(i)+1/6*(K1tg2+2*K2tg2+2*K3tg2+K4tg2);

tag(i+1)=tag(i)+1/6*(K1tag+2*K2tag+2*K3tag+K4tag);

tab(i+1)=tab(i)+1/6*(K1tab+2*K2tab+2*K3tab+K4tab);

tpot(i+1)=tpot(i)+1/6*(K1tpot+2*K2tpot+2*K3tpot+K4tpot);

tload(i+1)=tload(i)+1/6*(K1tload+2*K2tload+2*K3tload+K4tload);

end

%plot the solution

% plot(t,tf)

% hold on

% plot(t,tg1)

% hold on

% plot(t,tg2)

% hold on

% plot(t,tag)

% hold on

% plot(t,tpot)

% hold on

plot(t,tload)

Please help me to sort out this problem... Help will be appreciated.