MATLAB: Condition Not Working In Closed-Loop Simulation

I'm simulating the immune system response to pathogens and therapy. it's made so that if the level of organ damage reaches 50% or higher, the production of plasma cells stops. The condition works when I simulate it in open-loop but when it's in closed-loop, it doesn't.

The open-loop simulation, in which no therapy is applied, is run in this runimmune.m file:

% runimmune ? runs basic simulation of pathogen attack (no controls)
ipath=[1.5 2 2.57 3];
dT=.01;
tspan=0:dT:10;
u=zeros(size(tspan,2),4);
options=[];
handle=@f_of_x;
a11=1;a12=1;a23=1;a31=1;a42=1;b2=1;b3=1;
b1=-1;b4=-1;
a22=3;
a32=1.5;
a33=0.5;a41=0.5;
x2e=2;
x2e=2;x3e=a31*x2e/a32;
sty={'' '--' '-.' ':'};
for ii = 1:4
sol=ode45(handle,[0 10],[ipath(ii);2;2*x2e/3;0],options,dT,u);
subplot(2,2,1);plot(sol.x,sol.y(1,:),sty{ii},'linewidth',1.0);axis([0 10 0 6]);ylabel('pathogen conc.');grid on;hold on
subplot(2,2,2);plot(sol.x,sol.y(2,:),sty{ii},'linewidth',1.0);axis([0 10 0 6]);ylabel('plasma conc.');grid on;hold on
subplot(2,2,3);plot(sol.x,sol.y(3,:),sty{ii},'linewidth',1.0);axis([0 10 0 3]);ylabel('antibody conc.');grid on;hold on
subplot(2,2,4);refline(0,0.5);plot(sol.x,sol.y(4,:),sty{ii},'linewidth',1.0);axis([0 10 0 1]);ylabel('organ health');grid on;hold on
end
hold off
legend('Sub-Clinical','Clinical','Chronic','Lethal')

with f_of_x.m being

function dxdt=f_of_x(t,x,dT,u)
a11=1;a12=1;a23=1;a31=1;a42=1;b2=1;b3=1;
b1=-1;b4=-1;
a22=3;
a32=1.5;
a33=0.5;a41=0.5;
x2e=2;
n=floor(t/dT)+1;
dxdt=zeros(4,1);
dxdt(1)=(a11-a12*x(3))*x(1)+b1*u(n,1);
dxdt(2)=a21(x(4))*a22*x(1)*x(3)-a23*(x(2)-x2e)+b2*u(n,2);
dxdt(3)=a31*x(2)-(a32+a33*x(1))*x(3)+b3*u(n,3);
dxdt(4)=a41*x(1)-a42*x(4)+b4*u(n,4);
function y=a21(x)
if x>=0 && x<0.5
    y=cos(pi*x);
elseif x>=0.5
    y=0;
else
    error('negative values not allowed')
end

The condition is applied in the a21.m function (in the f_of_x.) The above code works fine as I wanted.

For the closed-loop, I have included the same condition but it's not called.

The closed-loop simulation is in runimmuneCL2.m

% runimmune CL  solves closed-loop non-linear optimal control for pathogen attack
clear; close all
ipath=[1.5 2 2.57 3]; % different initial pathogen concentrations (ICs)
sys_name='imresp2'; % define system name - could be handled more neatly
sys=imresp2; % get necessary system info
po=sys.f(:); % vectorize the F matrix
xe=sys.xe;   % this is the equilibrium value so we can re-construct actual concs.
To=0;Tf=10;  % set start and end times
options=[];  % options for ode solver (leaves at default)
dhandle=@dynamics;rhandle=@riccati; % set up handles for riccati eqn and dynamic eqn
thresh=1.e-3; % termination threshold
sty={'k' 'k--' 'k-.' 'k:'}; % set up line styles for plotting
% set up for computing values of control signal(s)
dT=.01;
tspan=To:dT:Tf;
old_u=zeros(1,length(tspan));
xo=[0 0 0 0]'; % set up IC - the first element will change depending on which level of pathogen attack is used
for ii = 1:length(ipath) % loop round the different pathogen attack scenarios
    term=100; % re-set termination initial value for each run
    solx=[];solp=[]; % erase previous values of P(t) and x(t)
    xo(1)=ipath(ii);solx=xo; % First iteration uses fixed value of xo.
    cnt=0; % counter
    while term>thresh
        solp=ode23s(rhandle,[0 Tf],po,options,sys_name,solx,Tf); % solve riccati ode
        solx=ode23s(dhandle,[0 Tf],xo,options,sys_name,solx,solp,Tf); % solve dynamics ode
        [u,J,Jf]=evalcontrol(tspan,sys_name,solx,solp,Tf); % compute control signal
        term=(norm(u-old_u))/Tf; % compute termination test value
        old_u=u; % over-write old controls with new
        disp(term); cnt=cnt+1;
    end
    disp(['No. Iterations = ' int2str(cnt)])
subplot(3,2,1);plot(solx.x,solx.y(1,:),sty{ii},'linewidth',1.0);axis([0 10 0 5]);ylabel('pathogen conc.');grid on;hold on
subplot(3,2,2);plot(solx.x,solx.y(2,:)+sys.xe(2),sty{ii},'linewidth',1.0);axis([0 10 0 10]);ylabel('plasma conc.');grid on;hold on
subplot(3,2,3);plot(solx.x,solx.y(3,:)+sys.xe(3),sty{ii},'linewidth',1.0);axis([0 10 0 5]);ylabel('antibody conc.');grid on;hold on
subplot(3,2,4);plot(solx.x,solx.y(4,:),sty{ii},'linewidth',1.0);axis([0 10 0 1]);ylabel('organ health');grid on;hold on
subplot(3,2,5);plot(tspan,u,sty{ii},'linewidth',1.0);grid on;hold on
subplot(3,2,6);plot(tspan,J,sty{ii},'linewidth',1.0);grid on;hold on
end

and imresp2.m is where the condition is written.

function sys=imresp2(t,x)
a11=1;a12=1;a23=1;a31=1;a42=1;b2=1;b3=1;
b1=-1;b4=-1;
a22=3;
a32=1.5;
a33=0.5;a41=0.5;
f11=1;f44=1;
q11=1;q44=1;
r11=1;
x2e=2;x3e=a31*x2e/a32;
if nargin==0;sys.f=[f11 0 0 0;0 0 0 0;0 0 0 0;0 0 0 f44];sys.xe=[0 x2e x3e 0]';
else
    sys.a=[((a11-a12*x3e)-a12*x(3)) 0 0 0;...
        a21(x(4))*a22*(x(3)+x3e) -a23 0 0;...
        -a33*x3e a31 -(a32+a33*x(1)) 0;...
        a41 0 0 -a42];
    sys.b=[b1 0 0 0]';
    sys.q=[q11 0 0 0;0 0 0 0;0 0 0 0;0 0 0 q44];
    sys.r=r11;
end
function y=a21(x)
if x>=0 && x<0.5
    y=cos(pi*x);
elseif x>=0.5
    y=0;
else
    error('negative values not allowed')

The exact same condition is written for a21. Can someone tell me why the condition is not applied in the closed-loop simulation?

(dynamics.m, riccati.m, and evalcontrol.m is included below.)

Thank you!

---

dynamics.m:

function dx=dynamics(t,x,sys_name,solx,solp,Tf)
if ~isstruct(solx);
    sys=feval(sys_name,t,solx); % first iteration sets to fixed value for all t, i.e. the IC.
else
    sys=feval(sys_name,t,deval(solx,Tf-t));
end
k=-sys.r\sys.b'*reshape(deval(solp,Tf-t),length(x),length(x));
dx=(sys.a+sys.b*k)*x;

riccati.m:

function dp=riccati(t,p,sys_name,solx,Tf)
n=sqrt(length(p));
P=reshape(p,n,n);
if ~isstruct(solx);
    sys=feval(sys_name,t,solx); % first iteration sets to fixed value for all t, i.e. the IC.
else
    sys=feval(sys_name,t,deval(solx,Tf-t));
end
dP=sys.q +P*sys.a +sys.a'*P-P*sys.b/sys.r*sys.b'*P;
dp=dP(:);

evalcontrol.m:

function [u,J,Jf]=evalcontrol(tspan,sys_name,solx,solp,Tf)
xf=deval(solx,Tf);
tmp=zeros(1,length(tspan));
u=zeros(1,length(tspan)); % sort out for >1 control input
for ii=1:length(tspan)
    t=tspan(ii);
    x=deval(solx,t);
    sys=feval(sys_name,t,x);
    u(:,ii)=-sys.r\sys.b'*reshape(deval(solp,Tf-t),length(x),length(x))*x;
    tmp(ii)=x'*sys.q*x+u(:,ii)'*sys.r*u(:,ii);
end
J=cumtrapz(tspan,tmp)/2;
sys=imresp2;
Jf=xf'*sys.f*xf/2;