function y = cap2_rubeola_back_forw
test = -1;
tf=3;
delta = 0.001;
M = 999;
t = linspace(0,tf,M+1);
h = tf/M;
h2 = h/2;
A(1)=10;
beta=36.5;
psi=30.417;
mu=0.65;
theta=0.65;
alpha=527.59;
k=0.136;
z=0.5;
delta(1)=0.7;
eta=0.134;
xi=0.12;
Q=10;
Q2=1;
tau=1.0;
x1=zeros(1,M+1);
x2=zeros(1,M+1);
x3=zeros(1,M+1);
x4=zeros(1,M+1);
u = zeros(1,M+1);
lambda1 = zeros(1,M+1);
lambda2 = zeros(1,M+1);
lambda3 = zeros(1,M+1);
lambda4 = zeros(1,M+1);
x1(1) = 0.0555;
x2(1) = 0.0003;
x3(1) = 0.0004;
x4(1) = 1;
while(test < 0)
oldu = u;
oldx1 = x1;
oldx2 = x2;
oldx3 = x3;
oldx4 = x4;
oldlambda1 = lambda1;
oldlambda2 = lambda2;
oldlambda3 = lambda3;
oldlambda4 = lambda4;
for i = 1:M
m11 = A(1)-beta*x1(i)*x2(i)-psi*x1(i)*x3(i)/(x3(i)+xi)-(mu+u(i))*x1(i)-theta*x4(i);
m12 = beta*x1(i)*x2(i)-psi*x1(i)*x3(i)/(x3(i)+xi)-(alpha+mu+k)*x2(i);
m13 = z*x2(i)-(delta+eta)*x3(i);
m14 = alpha*x2(i)-(theta+mu)*x4(i)+u(i)*x1(i);
m21 = A(1)-beta*(x1(i)+h2*m11)*(x2(i)+h2*m12)-psi*(x1(i)+h2*m11)*(x3(i)+h2*m13)/((x3(i)+h2*m13)+xi)-(mu+(0.5*(u(i) + u(i+1))))*(x1(i)+h2*m11)-theta*(x4(i)+h2*m14);
m22 = beta*(x1(i)+h2*m11)*(x2(i)+h2*m12)-psi*(x1(i)+h2*m11)*(x3(i)+h2*m13)/((x3(i)+h2*m13)+xi)-(alpha+mu+k)*(x2(i)+h2*m12);
m23 = z*(x2(i)+h2*m12)-(delta+eta)*(x3(i)+h2*m13);
m24 = alpha*(x2(i)+h2*m12)-(theta+mu)*(x4(i)+h2*m14)+((0.5*(u(i) + u(i+1))))*(x1(i)+h2*m11);
m31 = A(1)-beta*(x1(i)+h2*m21)*(x2(i)+h2*m22)-psi*(x1(i)+h2*m21)*(x3(i)+h2*m23)/((x3(i)+h2*m23)+xi)-(mu+(0.5*(u(i) + u(i+1))))*(x1(i)+h2*m21)-theta*(x4(i)+h2*m24);
m32 = beta*(x1(i)+h2*m21)*(x2(i)+h2*m22)-psi*(x1(i)+h2*m21)*(x3(i)+h2*m23)/((x3(i)+h2*m23)+xi)-(alpha+mu+k)*(x2(i)+h2*m22);
m33 = z*(x2(i)+h2*m22)-(delta+eta)*(x3(i)+h2*m23);
m34 = alpha*(x2(i)+h2*m22)-(theta+mu)*(x4(i)+h2*m24)+((0.5*(u(i) + u(i+1))))*(x1(i)+h2*m21);
m41 = A(1)-beta*(x1(i)+h2*m31)*(x2(i)+h2*m32)-psi*(x1(i)+h2*m31)*(x3(i)+h2*m33)/((x3(i)+h2*m33)+xi)-(mu+ u(i+1))*(x1(i)+h2*m31)-theta*(x4(i)+h2*m34);
m42 = beta*(x1(i)+h2*m31)*(x2(i)+h2*m32)-psi*(x1(i)+h2*m31)*(x3(i)+h2*m33)/((x3(i)+h2*m33)+xi)-(alpha+mu+k)*(x2(i)+h2*m32);
m43 = z*(x2(i)+h2*m32)-(delta+eta)*(x3(i)+h2*m33);
m44 = alpha*(x2(i)+h2*m32)-(theta+mu)*(x4(i)+h2*m34)+ u(i+1)*(x1(i)+h2*m31);
x1(i+1) = x1(i) + (h/6)*(m11 + 2*m21 + 2*m31 + m41);
x2(i+1) = x2(i) + (h/6)*(m12 + 2*m22 + 2*m32 + m42);
x3(i+1) = x3(i) + (h/6)*(m13 + 2*m23 + 2*m33 + m43);
x4(i+1) = x4(i) + (h/6)*(m14 + 2*m24 + 2*m34 + m44);
end
for i = 1:M
j = M + 2 – i;
n11 = Q – beta*x2*(lambda1(j)-lambda2(j))- psi*x3(j)/(x3(j)+ xi)*((lambda1(j)-lambda2(j)))-(mu + u(j));
n12 = Q2-lambda1(j)*beta*x1 + lambda2(j)*beta*x1-(alpha + mu + k)+lambda3(j)*z+lambda4(j)*alpha;
n13 = -lambda1(j)*psi*x3(j)*x1(j)/(x3(j)+ xi)^2;
n14 = -lambda4(j)*(theta+mu)+u(j)*x1(j);
n21 = Q – beta*(0.5*(x2(j)+x2(j-1)))*((lambda1(j)-h2*n11)-(lambda2(j)-h2*n12))- psi*(0.5*(x3(j)+x3(j-1)))/((x3(j)+x3(j-1))*xi)*((lambda1(j)-h2*n11)-(lambda2(j)-h2*n12))-(mu + u(j));
n22 = Q2-(lambda1(j)- h2*n11)*beta*(x1(j)+x1(j-1)) + (lambda2(j)-h2*n12)*beta*(x1(j)+x1(j-1))-(alpha + mu + k)+(lambda3(j)-h2*n13)*z+(lambda4(j)-h2*n14)*alpha;
n23 = -(lambda1(j)-h2*n11)*psi*(x3(j)+x3(j-1))*(x1(j)+x1(j-1))/(x3(j)+x3(j-1)+ xi)^2-(lambda2(j)-h2*n12)*psi*x3(j)*x1(j)/(x3(j)+ xi)^2-(lambda3(j)-h2*n13)*(delta+eta)+(lambda1(j)-h2*n11)*theta;
n24 = -(lambda4(j)-h2*n14)*(theta+mu)+u(j)*x1(j);
n31 = Q – beta*(0.5*(x2(j)+x2(j-1)))*((lambda1(j)-h2*n21)-(lambda2(j)-h2*n22))- psi*(0.5*(x3(j)+x3(j-1)))/((x3(j)+x3(j-1))*xi)*((lambda1(j)-h2*n21)-(lambda2(j)-h2*n22))-(mu + u(j));
n32 = Q2-(lambda1(j)- h2*n21)*beta*(x1(j)+x1(j-1)) + (lambda2(j)-h2*n22)*beta*(x1(j)+x1(j-1))-(alpha + mu + k)+(lambda3(j)-h2*n23)*z+(lambda4(j)-h2*n24)*alpha;
n33 = -(lambda1(j)-h2*n21)*psi*(x3(j)+x3(j-1))*(x1(j)+x1(j-1))/(x3(j)+x3(j-1)+ xi)^2-(lambda2(j)-h2*22)*psi*x3(j)*x1(j)/(x3(j)+ xi)^2-(lambda3(j)-h2*23)*(delta+eta)+(lambda1(j)-h2*n21)*theta;
n34 = -(lambda4(j)-h2*n24)*(theta+mu)+u(j)*x1(j);
n41 = Q – beta*(0.5*(x2(j)+x2(j-1)))*((lambda1(j)-h2*n31)-(lambda2(j)-h2*n32))- psi*(0.5*(x3(j)+x3(j-1)))/((x3(j)+x3(j-1))*xi)*((lambda1(j)-h2*n31)-(lambda2(j)-h2*n32))-(mu + u(j));
n42 = Q2-(lambda1(j)- h2*n31)*beta*(x1(j)+x1(j-1)) + (lambda2(j)-h2*n32)*beta*(x1(j)+x1(j-1))-(alpha + mu + k)+(lambda3(j)-h2*n33)*z+(lambda4(j)-h2*n34)*alpha;
n43 = -(lambda1(j)-h2*n31)*psi*(x3(j)+x3(j-1))*(x1(j)+x1(j-1))/(x3(j)+x3(j-1)+ xi)^2-(lambda2(j)-h2*32)*psi*x3(j)*x1(j)/(x3(j)+ xi)^2-(lambda3(j)-h2*33)*(delta+eta)+(lambda1(j)-h2*n31)*theta;
n44 = -(lambda4(j)-h2*n34)*(theta+mu)+u(j)*x1(j);
lambda1(j-1) = lambda1(j) – h/6*(n11 + 2*n21 + 2*n31 + n41);
lambda2(j-1) = lambda2(j) – h/6*(n12 + 2*n22 + 2*n32 + n42);
lambda3(j-1) = lambda3(j) – h/6*(n13 + 2*n23 + 2*n33 + n43);
lambda4(j-1) = lambda4(j) – h/6*(n14 + 2*n24 + 2*n34 + n44);
Best Answer