No output in my problem

Question

0 votes

I was wondering if somebody could help me about the following code which doesn't have output as a figure. Thanks in advance.

clc
clear
n = 1000;
a = 0;
b = 7;
h = (b-a)/n;
x(1) = 0.009;
y(1) = 0.001;
z(1)=0.99;
t(1)=a;
f=@(t,x1,x2,x3) -0.002*x1+0.018*x2-0.10*x1+(1-0.10)*(0.01+(0.99/(1+exp(-2*t+8)))+1)*((x3*x2)/(x1+x2+x3))+0.06*(1-0.05)*x3
g=@(t,x1,x2,x3) -0.018*x2+0.002*x1-0.10*x2+0.05*0.06*x3+0.10*(0.01+0.99/(1+exp(-2*t+8))+1)*((x3*x2)/(x1+x2+x3))
c=@(t,x1,x2,x3) -(0.01+0.99/(1+exp(-2*t+8))+1)*((x3*x2)/(x1+x2+x3))-0.06*x3+0.10*x2+0.10*x1
for i=1:n-1
    t(i+1) = a + i*h;
end
 disp('Step       t               x                      y              z');
 disp('----------------------------------------------------');
 fprintf('0       %2.2f       %2.16f       %2.16f      %2.16f\n', t(1),x(1),y(1),z(1));
 for k=1:n-1
f1 = f(t(k),x(k),y(k),z(k));
g1 = g(t(k),x(k),y(k),z(k));
c1 = c(t(k),x(k),y(k),z(k));
f2 = f(t(k)+(h/2), x(k)+(h/2)*f1 , y(k)+(h/2)*g1 ,z(k)+(h/2)*c1);
g2 = g(t(k)+(h/2), x(k)+(h/2)*f1 , y(k)+(h/2)*g1 ,z(k)+(h/2)*c1);
c2 = c(t(k)+(h/2), x(k)+(h/2)*f1 , y(k)+(h/2)*g1 ,z(k)+(h/2)*c1);
f3 = f(t(k)+(h/2), x(k)+(h/2)*f2 , y(k)+(h/2)*g2,z(k)+(h/2)*c2);
g3 = g(t(k)+(h/2), x(k)+(h/2)*f2 , y(k)+(h/2)*g2,z(k)+(h/2)*c2);
c3 = c(t(k)+(h/2), x(k)+(h/2)*f2 , y(k)+(h/2)*g2,z(k)+(h/2)*c2);
f4 = f(t(k)+h , x(k)+h*f3 , y(k)+h*g3,z(k)+h*c3);
g4 = g(t(k)+h , x(k)+h*f3 , y(k)+h*g3,z(k)+h*c3);
c4 = c(t(k)+h , x(k)+h*f3 , y(k)+h*g3,z(k)+h*c3);
x(k+1) = x(k) + (h/6)*(f1 + 2*f2 + 2*f3 + f4);
y(k+1) = y(k) + (h/6)*(g1+ 2*g2 + 2*g3 + g4);
z(k+1) = z(k) + (h/6)*(c1+ 2*c2 + 2*c3 + c4);
fprintf('%2f       %2.2f       %2.16f       %2.16f     %2.16f\n', k, t(k+1),x(k+1),y(k+1), z(k+1));
 end
 
 n = 1000;
a = 0;
b = 7;
h = (b-a)/n;
p(1) = 0;
q(1) = 0;
r(1) = 0;
t(1) = b;
f=@(t,x1,x2,x3) 0.018*(x1-x2)+0.10*(x1-x3)+((x3-0.10*x1-(1-0.10)*x2)*(0.01+(0.99)/(1+exp(-2*t+8))+1)*z(k)*(x(k)+z(k)))/((x(k)+y(k)+z(k))^2)
g=@(t,x1,x2,x3) 0.002*(x2-x1)+0.10*(x2-x3)-((x3-0.10*x1-(1-0.10)*x2)*(0.01+(0.99)/(1+exp(-2*t+8))+1)*z(k)*y(k))/((x(k)+y(k)+z(k))^2)
c=@(t,x1,x2,x3) -1+(x3-0.05*x1-(1-0.05)*x2)*0.06+((x3-0.10*x1-(1-0.10)*x2)*(0.01+(0.99)/(1+exp(-2*t+8))+1)*y(k)*(x(k)+y(k)))/((x(k)+y(k)+z(k))^2)
for i=1:n-1
    t(i+1) = b - i*h;
end
 disp('Step       t               x                      y              z');
 disp('----------------------------------------------------');
 fprintf('0       %2.2f       %2.16f       %2.16f      %2.16f\n', t(1),p(1),q(1),r(1));
 for k=1:n-1
f1 = f(t(k),p(k),q(k),r(k));
g1 = g(t(k),p(k),q(k),r(k));
c1 = c(t(k),p(k),q(k),r(k));
f2 = f(t(k)-(h/2), p(k)-(h/2)*f1 , q(k)-(h/2)*g1 ,r(k)-(h/2)*c1);
g2 = g(t(k)-(h/2), p(k)-(h/2)*f1 , q(k)-(h/2)*g1 ,r(k)-(h/2)*c1);
c2 = c(t(k)-(h/2), p(k)-(h/2)*f1 , q(k)-(h/2)*g1 ,r(k)-(h/2)*c1);
f3 = f(t(k)-(h/2), p(k)-(h/2)*f2 , q(k)-(h/2)*g2, r(k)-(h/2)*c2);
g3 = g(t(k)-(h/2), p(k)-(h/2)*f2 , q(k)-(h/2)*g2, r(k)-(h/2)*c2);
c3 = c(t(k)-(h/2), p(k)-(h/2)*f2 , q(k)-(h/2)*g2, r(k)-(h/2)*c2);
f4 = f(t(k)-h , p(k)-h*f3 , q(k)-h*g3,r(k)-h*c3);
g4 = g(t(k)-h , p(k)-h*f3 , q(k)-h*g3,r(k)-h*c3);
c4 = c(t(k)-h , p(k)-h*f3 , q(k)-h*g3,r(k)-h*c3);
p(k+1) = p(k) - (h/6)*(f1 + 2*f2 + 2*f3 + f4);
q(k+1) = q(k) - (h/6)*(g1+ 2*g2 + 2*g3 + g4);
r(k+1) = r(k) - (h/6)*(c1+ 2*c2 + 2*c3 + c4);
fprintf('%2f       %2.2f       %2.16f       %2.16f     %2.16f\n', k, t(k+1),p(k+1),q(k+1), r(k+1));
 end
 
 figure
plot(t,min(max(0,((r-(1/10)*p-(1-(1/10))*q).*z.*y)/(2*(1/100)*(x+y+z))),1))
axis([0 7 0 1]);
xlabel('time')
ylabel('Control U2')