Codding for Numerical Analysis using control

function y = code9(shy,beta,alpha,epsilon,gama,mu,rI,rH,c,d,S0,Q0,I0,H0,R0,,u1,u2,u3,A,B,C,T)
shy=.0000000032;
beta=.03;
alpha=.1;
epsilon=.01;
d=.005;
gama=.6;
rI=.001;
c=0.4;
mu=.1;
rH=.7;

S0 = 10000;
Q0 = 6000;
I0 = 70;
H0 = 30;
R0 = 20;
A = 60;
B = 80;
C=20
T = 10;
test = -1;
delta = 0.001;
N = 1000;
t=linspace(0,T,N+1);
h=T/N;
h2 = h/2;

S=zeros(1,N+1);
S(1)=S0;
lamS=zeros(1,N+1);

Q=zeros(1,N+1);
Q(1)=Q0;
lamQ=zeros(1,N+1);

I=zeros(1,N+1);
I(1)=I0;
lamI=zeros(1,N+1);

H=zeros(1,N+1);
H(1)=H0;
lamH=zeros(1,N+1);

R=zeros(1,N+1);
R(1)=R0;
lamR=zeros(1,N+1);
u1=zeros(1,N+1);
u2=zeros(1,N+1);
u3=zeros(1,N+1);

while(test < 0)

oldu1 = u1;
oldu2 = u2;
oldu3 = u3;
oldS = S;
oldQ = Q;
oldI = I;
oldH = H;
oldR = R;

oldlamS = lamS;
oldlamQ = lamQ;
oldlamI = lamI;
oldlamH = lamH;
oldlamR = lamR;

for i=1:N
k11 = shy-beta*S(i)*I(i)-gama*S(i)*Q(i)—u1(i)*S(i)+(1-alpha)*u2(i)*Q(i)+epsilon*R(i)-d*S(i);

k21 = beta*S(i)*I(i)+ gama*S(i)*Q(i)-alpha*u2(i)*Q(i)-( 1-alpha)*u2(i)*Q(i)-d*Q(i);

k31 = alpha*u2(i)*Q(i)-rI*I(i)-c*I(i)-mu*I(i)-d*I(i);

k41 = c*I(i)-(rH+mu+d)*H(i)-u3(i)*H(i);

k51 = rH*H(i)+rI*I(i)-epsilon*R(i)-d*R(i)+u1(i)*S(i)+u3(i)*H(i);
k12=shy-beta*(S(i)+h2*k11)*(I(i)+h2*k31)-gama*(S(i)+h2*k11)*(Q(i)+h2*k21)-.5*(u1(i)
+u1(i+1))*(S(i)+h2*k11)+(1-alpha)*(0.5*(u2(i)+u2(i+1))*(Q(i)+h2*k21)+epsilon*(R(i)+h2*k51)-d*(S(i)
+h2*k11);

k22=beta*(S(i)+h2*k11)*(I(i)+h2*k31)+gama*(S(i)+h2*k11)*(Q(i)+h2*k21)-alpha*(0.5*(u2(i)
+u2(i+1))*(Q(i)+h2*k21)-(1-alpha)*(0.5*(u2(i)+u2(i+1))*(Q(i)+h2*k21)-d*(Q(i)+h2*k21);
k32=alpha*0.5(u2(i)+u2(i+1))*(Q(i)+h2*k21)-rI*(I(i)+h2*k31)-c*(I(i)+h2*k31)-mu*(I(i)+h2*k31)-
d*(I(i)+h2*k31);
k42= c*(I(i)+h2*k31)-(rH+mu+d)*(H(i)+h2*k41)-0.5*(u3(i)+u3(i+1))*(H(i)+h2*k41);

k52=rH*(H(i)+h2*k41)+rI*(I(i)+h2*k31)-epsilon*(R(i)+h2*k51)-d*(R(i)+h2*k51)+0.5*(u1(i)
+u1(i+1))*(S(i)+h2*k11)+0.5*(u3(i)+u3(i+1))*(H(i)+h2*k41);
k13=shy-beta*(S(i)+h2*k12)*(I(i)+h2*k32)-gama*(S(i)+h2*k12)*(Q(i)+h2*k22)-.5*(u1(i)
+u1(i+1))*(S(i)+h2*k12)+(1-alpha)*(0.5*(u2(i)+u2(i+1))*(Q(i)+h2*k22)+epsilon*(R(i)+h2*k52)-d*(S(i)
+h2*k12);
 k23=beta*(S(i)+h2*k12)*(I(i)+h2*k32)+gama*(S(i)+h2*k12)*(Q(i)+h2*k22)-alpha*(0.5*(u2(i)
+u2(i+1))*(Q(i)+h2*k22)-(1-alpha)*(0.5*(u2(i)+u2(i+1))*(Q(i)+h2*k22)-d*(Q(i)+h2*k22);

k33=alpha*0.5(u2(i)+u2(i+1))*(Q(i)+h2*k22)-rI*(I(i)+h2*k32)-c*(I(i)+h2*k32)-mu*(I(i)+h2*k32)-
d*(I(i)+h2*k32);

k43= c*(I(i)+h2*k32)-(rH+mu+d)*(H(i)+h2*k42)-0.5*(u3(i)+u3(i+1))*(H(i)+h2*k42);

k53=rH*(H(i)+h2*k42)+rI*(I(i)+h2*k32)-epsilon*(R(i)+h2*k52)-d*(R(i)+h2*k52)+0.5*(u1(i)
+u1(i+1))*(S(i)+h2*k12)+0.5*(u3(i)+u3(i+1))*(H(i)+h2*k42);
k14=shy-beta*(S(i)+h2*k13)*(I(i)+h2*k33)-gama*(S(i)+h2*k13)*(Q(i)+h2*k23)-u1(i+1))*(S(i)
+h2*k13)+(1-alpha)*u2(i+1)*(Q(i)+h2*k23)+epsilon*(R(i)+h2*k53)-d*(S(i)+h2*k13);

k24=beta*(S(i)+h2*k13)*(I(i)+h2*k33)+gama*(S(i)+h2*k13)*(Q(i)+h2*k23)-alpha*(u2(i+1))*(Q(i)
+h2*k23)-(1- alpha)*(u2(i+1))*(Q(i)+h2*k23)-d*(Q(i)+h2*k23);

k34=alpha*u2(i+1)*(Q(i)+h2*k23)-rI*(I(i)+h2*k33)-c*(I(i)+h2*k33)-mu*(I(i)+h2*k33)-d*(I(i)
+h2*k33);

k44= c*(I(i)+h2*k33)-(rH+mu+d)*(H(i)+h2*k43)-u3(i+1)*(H(i)+h2*k43);

k54=rH*(H(i)+h2*k43)+rI*(I(i)+h2*k33)-epsilon*(R(i)+h2*k53)-d*(R(i)+h2*k53)+u1(i+1)*(S(i)
+h2*k13)+u3(i+1)*(H(i)+h2*k43);

S(i+1) = S(i) + (h/6)*(k11 + 2*k12 + 2*k13 + k14);

Q(i+1) = Q(i) + (h/6)*(k21 + 2*k22 + 2*k23 + k24);
I(i+1) = I(i) + (h/6)*(k31 + 2*k32 + 2*k33 + k34);
H(i+1) = H(i) + (h/6)*(k41 + 2*k42 + 2*k43 + k44);
R(i+1) = R(i) + (h/6)*(k51 + 2*k52 + 2*k53 + k54);

end

for i = 1:N
j = N + 2 - i;
k11 = (beta*I(j)+u1(j)+d+gama*Q(j))*lamS(j)-(beta*I(j)-gama*Q(j))*lamQ(j)-
u1(j)*lamR(j);
k21=(-(1-alpha)*u2(j)+gama*S(j))*lamS(j)-(-alpha+gama*S(j)-(1-
alpha)*u2(j)-d)*lamQ(j)-alpha*u2(j)*lamI(j);
k31=1+beta*S(j)*lamS(j)-beta*S(j)*lamQ(j)+(rI+c+mu+d)*lamI(j)-
c*lamH(j)-rI*lamR(j);
k41=(rH+mu+d+u3(j))*lamH(j)-(rH+u3(j))*lamR(j);
k51=-epsilon*lamS(j)+(epsilon+d)*lamR(j);
k12=(beta*0.5(I(j)+I(j-1))+0.5*(u1(j)+u1(j-1))+d+gama*0.5(Q(j)+Q(j-
1)))*(lamS(j)-h2*k11)-(beta*0.5*(I(j)+I(j-1))-gama*0.5*(Q(j)+Q(j-1))*(lamQ(j)-
h2*k21)-0.5*(u1(j)+u1(j-1))*(lamR(j)-h2*k51);
k22 = (-(1-alpha)*0.5*(u2(j)+u2(j-1))+gama*0.5*(S(j)+S(j-1))*(lamS(j)-
h2*k11)-(-alpha+gama*0.5*(S(j)+S(j-1))-(1-alpha)*0.5*(u2(j)+u2(j-1))-
d)*(lamQ(j)-h2*k21)-alpha*0.5*(u2(j)+u2(j-1))*(lamI(j)-h2*k31);
 k32=1+beta*0.5*(S(j)+S(j-1))*(lamS(j)-h2*k11)-beta*0.5*(S(j)+S(j-
1))*(lamQ(j)-h2*k21)+(rI+c+mu+d)*(lamI(j)-h2*k31)-c*(lamH(j)-h2*k41)-
rI*(lamR(j)-h2*k51);
k42=(rH+mu+d+0.5*(u3(j)+u3(j-1)))*(lamH(j)-h2*k41)-(rH+0.5*(u3(j)
+u3(j-1)))*(lamR(j)-h2*k51);
k52=- epsilon*(lamS(j)-h2*k11)+(epsilon+d)*(lamR(j)-h2*k51);
k13==(beta*0.5(I(j)+I(j-1))+0.5*(u1(j)+u1(j-1))+d+gama*0.5(Q(j)+Q(j-
1)))*(lamS(j)-h2*k12)-(beta*0.5*(I(j)+I(j-1))-gama*0.5*(Q(j)+Q(j-1))*(lamQ(j)-
h2*k22)-0.5*(u1(j)+u1(j-1))*(lamR(j)-h2*k52);

k23 =(-(1-alpha)*0.5*(u2(j)+u2(j-1))+gama*0.5*(S(j)+S(j-1))*(lamS(j)-
h2*k12)-(-alpha+gama*0.5*(S(j)+S(j-1))-(1-alpha)*0.5*(u2(j)+u2(j-1))-
d)*(lamQ(j)-h2*k22)-alpha*0.5*(u2(j)+u2(j-1))*(lamI(j)-h2*k32);
k33=1+beta*0.5*(S(j)+S(j-1))*(lamS(j)-h2*k12)-beta*0.5*(S(j)+S(j-
1))*(lamQ(j)-h2*k22)+(rI+c+mu+d)*(lamI(j)-h2*k32)-c*(lamH(j)-h2*k42)-
rI*(lamR(j)-h2*k52);

k43=(rH+mu+d+0.5*(u3(j)+u3(j-1)))*(lamH(j)-h2*k42)-(rH+0.5*(u3(j)
+u3(j-1)))*(lamR(j)-h2*k52);

k53=- epsilon*(lamS(j)-h2*k12)+(epsilon+d)*(lamR(j)-h2*k52);

k14==(beta*I(j-1)+u1(j-1)+d+gama*Q(j-1))*(lamS(j)-h2*k13)-(beta*I(j-
1)-gama*Q(j-1)*(lamQ(j)-h2*k23)-u1(j-1)*(lamR(j)-h2*k5);

k24 = (-(1-alpha)*u2(j-1)+gama*S(j-1)*(lamS(j)-h2*k13)-(-
alpha+gama*S(j-1)-(1-alpha)*u2(j-1)-d)*(lamQ(j)-h2*k23)-alpha*u2(j-
1)*(lamI(j)-h2*k33);
k34=1+beta*S(j-1)*(lamS(j)-h2*k13)-beta*S(j-1)*(lamQ(j)-h2*k23)+
(rI+c+mu+d)*(lamI(j)-h2*k33)-c*(lamH(j)-h2*k43)-rI*(lamR(j)-h2*k53);
k44=(rH+mu+d+u3(j-1))*(lamH(j)-h2*k43)-(rH+u3(j-1))*(lamR(j)-h2*k53);

k54=- epsilon*(lamS(j)-h2*k13)+(epsilon+d)*(lamR(j)-h2*k53);

lamS(j-1) = lamS(j) - (h/6)*(k11 + 2*k12 + 2*k13 + k14);

lamQ(j-1) = lamQ(j) - (h/6)*(k21 + 2*k22 + 2*k23 + k24);
lamI(j-1) = lamI(j) - (h/6)*(k31 + 2*k32 + 2*k33 + k34);
lamH(j-1) = lamH(j) - (h/6)*(k41 + 2*k42 + 2*k43 + k44);
lamR(j-1) = lamR(j) - (h/6)*(k51 + 2*k52 + 2*k53 + k54);

end

tempp = (lamS-lamR)./(2*A);
u11 = max(0, min(1,tempp));
u1 = (u11 + oldu1)/2;
tempf = (alpha*Q(lamS-lamI)+Q(lamQ-lamS))/(2*beta)
u21 = max(0, min(1,tempf));
 u2 = (u21 + oldu2)/2;
tempq = (H*(lamH-lamR).)/(2*C);
u31 = max(0, min(1,tempq));
u3 = (u31 + oldu3)/2;

temp1 = delta*sum(abs(u1)) - sum(abs(oldu1- u1));

temp2 = delta*sum(abs(u2)) - sum(abs(oldu2 – u2));
temp3 = delta*sum(abs(u3)) - sum(abs(oldu3 – u3));
temp4 = delta*sum(abs(S)) - sum(abs(oldS - S));
temp5 = delta*sum(abs(Q)) - sum(abs(oldQ - Q));
temp6 = delta*sum(abs(I)) - sum(abs(oldI - I));
temp7 = delta*sum(abs(H)) - sum(abs(oldH - H));
temp8 = delta*sum(abs(R)) - sum(abs(oldR - R));
temp9 = delta*sum(abs(lamS)) - sum(abs(oldlamS - lamS));
temp10 = delta*sum(abs(lamQ)) - sum(abs(oldlamQ - lamQ));
temp11 = delta*sum(abs(lamI)) - sum(abs(oldlamI - lamI));
temp12 = delta*sum(abs(lamH)) - sum(abs(oldlamH - lamH));
temp13 = delta*sum(abs(lamR)) - sum(abs(oldlamR - lamR));

test = min(temp1, min(temp2, min(temp3, min(temp4, min(temp5, min(temp6,

min(temp7,min(temp8,min(temp9,min(temp10,min(temp11,min(temp12,temp13))))))))))));
end

y(1,:) = t;
y(2,:) = S;
y(3,:) = Q;
y(4,:) = I;
y(5,:) = H;
y(6,:) = R;
y(7,:) = up;
y(8,:) = uf;
subplot(3,2,1);plot(y(1,:),y(2,:),'c-','LineWidth',1.2)
subplot(3,2,1);xlabel('Time(Days)')
subplot(3,2,1);ylabel('Susceptible(S)')
subplot(3,2,2);plot(y(1,:),y(3,:),'m','LineWidth',1.2)
subplot(3,2,2);xlabel('Time(Days)')
subplot(3,2,2);ylabel('Quaraintained (Q)')
subplot(3,2,3);plot(y(1,:),y(4,:),'g-','LineWidth',1.2)
subplot(3,2,3);xlabel('Time(Days)')
subplot(3,2,3);ylabel('Infected (I)')
subplot(3,2,4);plot(y(1,:),y(5,:),'y','LineWidth',2)
subplot(3,2,4);xlabel('Time(Days)')
subplot(3,2,4);ylabel('Hospitalized (Is)')
subplot(3,2,5);plot(y(1,:),y(6,:),'b','LineWidth',1.2)
subplot(3,2,5);xlabel('Time(Days)')
subplot(3,2,5);ylabel('Recovered (R)')
 Codding for Numerical Analysis without control

clear all;clc;

%Number of Iterations 'iterations' multiplication of 6

iterations=1152;
loop=100;

% value of parameters
shy=32000;beta=.011;alpha=.5;epsylon=.04;d=.5;gama=.0000005;rI=.1;c=0.32;m=.99
;rH=.6;

%make the value of h

% upper limit 'b1' and lower limit 'a1'
up1=20;
low1=0;
h=(up1-low1)/iterations;
fS=@(t,S,Q,I,H,R) shy-beta*S*I-gama*S*Q-d*S+epsylon*R+(1-alpha)*Q;
fQ=@(t,S,Q,I,H,R) gama*S*Q+beta*S*I-alpha*Q-(1-alpha)*Q-d*Q;
fI=@(t,S,Q,I,H,R) alpha*Q-rI*I-c*I-m*I-d*I;
fH=@(t,S,Q,I,H,R) c*I-rH*H-m*H-d*H;
fR=@(t,S,Q,I,H,R) rI*I+rH*H-epsylon*R-d*R;

for i1=1:loop
% initial value
t(1)=0;
S(1)=10000;% at t=0
Q(1)=6000;% at t=0
I(1)=3000;% at t=0
H(1)=200;% at t=0
R(1)=20;% at t=0

%code using runge kutta method for state

for i=1:iterations
k11=fS(t(i),S(i),Q(i),I(i),H(i),R(i));
k21=fQ(t(i),S(i),Q(i),I(i),H(i),R(i));
k31=fI(t(i),S(i),Q(i),I(i),H(i),R(i));
k41=fH(t(i),S(i),Q(i),I(i),H(i),R(i));
k51=fR(t(i),S(i),Q(i),I(i),H(i),R(i));

k12=fS(t(i)+h/2,S(i)+h*k11/2,Q(i)+h*k21/2,I(i)+h*k31/2,H(i)
+h*k41/2,R(i)+h*k51/2);
k22=fQ(t(i)+h/2,S(i)+h*k11/2,Q(i)+h*k21/2,I(i)+h*k31/2,H(i)
+h*k41/2,R(i)+h*k51/2);
k32=fI(t(i)+h/2,I(i)+h*k11/2,Q(i)+h*k21/2,I(i)+h*k31/2,H(i)
+h*k41/2,R(i)+h*k51/2);
k42=fH(t(i)+h/2,S(i)+h*k11/2,Q(i)+h*k21/2,I(i)
+h*k31/2,H(i)+h*k41/2,R(i)+h*k51/2);
k52=fR(t(i)+h/2,S(i)+h*k11/2,Q(i)+h*k21/2,I(i)
+h*k31/2,H(i)+h*k41/2,R(i)+h*k51/2);
 k13=fS(t(i)+h/2,S(i)+h*k12/2,Q(i)+h*k22/2,I(i)+h*k32/2,H(i)
+h*k42/2,R(i)+h*k52/2);
k23=fQ(t(i)+h/2,S(i)+h*k12/2,Q(i)+h*k22/2,I(i)+h*k32/2,H(i)
+h*k42/2,R(i)+h*k52/2);
k33=fI(t(i)+h/2,S(i)+h*k12/2,Q(i)+h*k22/2,I(i)+h*k32/2,H(i)
+h*k42/2,R(i)+h*k52/2);
k43=fH(t(i)+h/2,S(i)+h*k12/2,Q(i)+h*k22/2,I(i)
+h*k32/2,H(i)+h*k42/2,R(i)+h*k52/2);
k53=fR(t(i)+h/2,S(i)+h*k12/2,Q(i)+h*k22/2,I(i)
+h*k32/2,H(i)+h*k42/2,R(i)+h*k52/2);

k14=fS(t(i)+h/2,S(i)+h*k13/2,Q(i)+h*k23/2,I(i)+h*k33/2,H(i)
+h*k43/2,R(i)+h*k53/2);
k24=fQ(t(i)+h/2,S(i)+h*k13/2,Q(i)+h*k23/2,I(i)+h*k33/2,H(i)
+h*k43/2,R(i)+h*k53/2);
k34=fI(t(i)+h/2,S(i)+h*k13/2,Q(i)+h*k23/2,I(i)+h*k33/2,H(i)
+h*k43/2,R(i)+h*k53/2);
k44=fH(t(i)+h/2,S(i)+h*k13/2,Q(i)+h*k23/2,I(i)
+h*k33/2,H(i)+h*k43/2,R(i)+h*k53/2);
k54=fR(t(i)+h/2,S(i)+h*k13/2,Q(i)+h*k23/2,I(i)
+h*k33/2,H(i)+h*k43/2,R(i)+h*k53/2);

delS1=h*(k11+2*k12+2*k13+k14)/6;
delQ1=h*(k21+2*k22+2*k23+k24)/6;
delI1=h*(k31+2*k32+2*k33+k34)/6;
delH1=h*(k41+2*k42+2*k43+k44)/6;
delR1=h*(k51+2*k52+2*k53+k54)/6;

t(i+1)=t(i)+h;
S(i+1)=S(i)+delS1;
Q(i+1)=Q(i)+delQ1;
I(i+1)=I(i)+delI1;
H(i+1)=H(i)+delH1;
R(i+1)=R(i)+delR1;

end

end
%

%
S(iterations+1)
Q(iterations+1)
I(iterations+1)
H(iterations+1)
R(iterations+1)
%

%time spacing for graph display function, control by "con1"

con1=h;
t1=0;
for i=1:iterations+1
t(i)=t1;
t1=t1+con1;
end
%t=[0,0.05,.1,.15,.2,.25,.3,.35,.4,.45,.5,.55,.6,.65,.7,.75,.8,.85,.9,.95,1];

figure
plot(t,S,'b',t,Q,'y',t,I,'r',t,H,'m',t,R,'g','linewidth',3);
xlabel('Time (days)','FontSize',10), ylabel('Population ');