I have a code for optimal control of a deterministic SEIR model. I am unable to code the stochastic version of the model.

test  = -1;         % Test variable; as long as variable is negative ...the while loops keeps repeating
t0    = 0;          % Initial time
tf    = 100;        % Final time
Delta = 0.00001;    % Accepted tollerance
M     = 999;        % Number of nodes
t     = linspace(t0, tf, M+1);  % Time variable where linspace creates M+1 ...equally spaced nodes between t0 and tf, including t0 and tf.
h     = tf/M;       % Spacing between nodes
h2    = h/2;        % Spacing equal to 2 for Runge-Kutta method
% MODEL PARAMETERS
alpha = 1.5;
mu = 0.05;
b1 = 0.05;
b2 = 0.03;
lambda = 1.2;
epsilon2 = 0.55;
epsilon1 = 0.01;
gamma = 0.2;
% WEIGHT FACTORS
k1 = 0.05;             % Weight on threatened population
k2 = 0.05;             % Weigth on deceased population
l1 = 0.05;
l2 = 0.05;
m1 = 0.025;
m2 = 0.025;
m3 = 0.025;
m4 = 0.025;
% INITIAL CONDITIONS MODEL
S  = zeros(1, M+1);  % Susceptible
E  = zeros(1, M+1);  % Exposed
I  = zeros(1, M+1);  % Infected
R  = zeros(1, M+1);  % Recovered
S(1) = 99999;
E(1) = 100;
I(1) = 100;
R(1) = 50;
% INITIAL GUESS FOR OPTIMAL CONTROL INPUT
u1 = zeros(1, M+1);     % Control input for government intervention
u2 = zeros(1, M+1);     % Control input for testing
% INITIAL CONDITIONS ADOINT SYSTEM
lambda1 = zeros(1, M+1);
lambda2 = zeros(1, M+1);
lambda3 = zeros(1, M+1);
lambda4 = zeros(1, M+1);
lambda1(M+1) = 0;
lambda2(M+1) = 0;
lambda3(M+1) = 0;
lambda4(M+1) = 0;
% FORWARD-BACKWARD SWEEP METHOD
loopcnt = 0; % Count number of loops
while(test < 0)
    loopcnt = loopcnt + 1;
    oldu1 = u1;
    oldu2 = u2;
    oldS = S;
    oldE = E;
    oldI = I;
    oldR = R;
    oldlambda1 = lambda1;
    oldlambda2 = lambda2;
    oldlambda3 = lambda3;
    oldlambda4 = lambda4;
    
    % SYSTEM DYNAMICCS
    for i = 1:M
        
        if i > round(tau1 / h)
            E_tau1 = E(i-round(tau1/h));
        else
            E_tau1 = E(1);
        end
        
        if i > round(tau2 / h)
            I_tau2 = I(i-round(tau2/h));
        else
            I_tau2 = I(1);
        end
        m11 = alpha - b1 * (1-u1(i)) * S(i) * E(i) / (S(i)+E(i)+I(i)+R(i)) - b2 * (1-u1(i)) * S(i) * I(i) / (S(i)+E(i)+I(i)+R(i))- mu * S(i);
        m12 = b1 * (1-u1(i)) * S(i) * E(i) / (S(i)+E(i)+I(i)+R(i)) + b2 * (1-u1(i)) * S(i) * E(i) / (S(i)+E(i)+I(i)+R(i)) - lambda * E_tau1 - (epsilon1 + mu) * E(i);
        m13 =  lambda * E_tau1 - epsilon2 * (1+u2(i)) * I_tau2 - (gamma + mu) * I(i);
        m14 = epsilon1 * (1+u2(i)) * E(i) + epsilon2 * I_tau2 - mu * R(i);
        
        m21 = alpha - b1 * ( 1 - 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m11) * (E(i)+h2*m12) / ((S(i)+h2*m11)+(E(i)+h2*m12)+(I(i)+h2*m13)+(R(i)+h2*m14)) - b2 * ( 1 - 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m11) * (I(i)+h2*m13) / ((S(i)+h2*m11)+(E(i)+h2*m12)+(I(i)+h2*m13)+(R(i)+h2*m14)) - mu * (S(i)+h2*m11);
        m22 = b1 * ( 1 - 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m11) * (E(i)+h2*m12) / ((S(i)+h2*m11)+(E(i)+h2*m12)+(I(i)+h2*m13)+(R(i)+h2*m14)) + b2 * ( 1 - 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m11) * (I(i)+h2*m13) / ((S(i)+h2*m11)+(E(i)+h2*m12)+(I(i)+h2*m13)+(R(i)+h2*m14)) - lambda * E_tau1 - (epsilon1 + mu) * (E(i)+h2*m12);
        m23 =  lambda * E_tau1 - epsilon2 * ( 1 + 0.5*(u2(i)+u2(i+1))) * I_tau2 - (gamma + mu) * (I(i)+h2*m13);
        m24 = epsilon1 * ( 1 + 0.5*(u2(i)+u2(i+1))) * (E(i)+h2*m12) + epsilon2 * I_tau2 - mu * (R(i)+h2*m14);
        m31 = alpha - b1 * ( 1 - 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m21) * (E(i)+h2*m22) / ((S(i)+h2*m21)+(E(i)+h2*m22)+(I(i)+h2*m23)+(R(i)+h2*m24)) - b2 * ( 1 - 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m21) * (I(i)+h2*m23) / ((S(i)+h2*m21)+(E(i)+h2*m22)+(I(i)+h2*m23)+(R(i)+h2*m24)) - mu * (S(i)+h2*m21);
        m32 = b1 * ( 1 - 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m21) * (E(i)+h2*m22) / ((S(i)+h2*m21)+(E(i)+h2*m22)+(I(i)+h2*m23)+(R(i)+h2*m24)) + b2 * ( 1 - 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m21) * (I(i)+h2*m23) / ((S(i)+h2*m21)+(E(i)+h2*m22)+(I(i)+h2*m23)+(R(i)+h2*m24)) - lambda * E_tau1 - (epsilon1 + mu) * (E(i)+h2*m22);
        m33 =  lambda * E_tau1 - epsilon2 * ( 1 + 0.5*(u2(i)+u2(i+1))) * I_tau2 - (gamma + mu) * (I(i)+h2*m23);
        m34 = epsilon1 * ( 1 + 0.5*(u2(i)+u2(i+1))) * (E(i)+h2*m22) + epsilon2 * I_tau2 - mu * (R(i)+h2*m24);
        m41 = alpha - b1 * ( 1 - 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m31) * (E(i)+h2*m32) / ((S(i)+h2*m31)+(E(i)+h2*m32)+(I(i)+h2*m33)+(R(i)+h2*m34)) - b2 * ( 1 - 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m31) * (I(i)+h2*m33) / ((S(i)+h2*m31)+(E(i)+h2*m32)+(I(i)+h2*m33)+(R(i)+h2*m34)) - mu * (S(i)+h2*m31);
        m42 = b1 * ( 1 - 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m31) * (E(i)+h2*m32) / ((S(i)+h2*m31)+(E(i)+h2*m32)+(I(i)+h2*m33)+(R(i)+h2*m34)) + b2 * ( 1 - 0.5*(u1(i)+u1(i+1))) * (S(i)+h2*m31) * (I(i)+h2*m33) / ((S(i)+h2*m31)+(E(i)+h2*m32)+(I(i)+h2*m33)+(R(i)+h2*m34)) - lambda * E_tau1 - (epsilon1 + mu) * (E(i)+h2*m32);
        m43 =  lambda * E_tau1 - epsilon2 * ( 1 + 0.5*(u2(i)+u2(i+1))) * I_tau2 - (gamma + mu) * (I(i)+h2*m33);
        m44 = epsilon1 * ( 1 + 0.5*(u2(i)+u2(i+1))) * (E(i)+h2*m32) + epsilon2 * I_tau2 - mu * (R(i)+h2*m34);
        
        
        S(i+1) = S(i) + (h/6)*(m11 + 2*m21 + 2*m31 + m41);
        E(i+1) = E(i) + (h/6)*(m12 + 2*m22 + 2*m32 + m42);
        I(i+1) = I(i) + (h/6)*(m13 + 2*m23 + 2*m33 + m43);
        R(i+1) = R(i) + (h/6)*(m14 + 2*m24 + 2*m34 + m44);
    end
    
    % ADJOINT SYSTEM
    for i = 1:M % From initial to final value
        j       = M + 2 - i;    % From final value to initial value
        n11 = -m1 * S(j) + lambda1(j) * b1 * (1-u1(j)) * E(j) / (S(j)+E(j)+I(j)+R(j)) + lambda1(j) * b2 * (1-u1(j)) * I(j) / (S(j)+E(j)+I(j)+R(j)) + mu * lambda1(j) - b1 * lambda2(j) * (1-u1(j)) * E(j) / (S(j)+E(j)+I(j)+R(j)) - b2 * lambda2(j) * (1-u1(j)) * I(j) / (S(j)+E(j)+I(j)+R(j));
        n12 = -k1 - m2 * E(j) + b1 * lambda1(j) * (1-u1(j)) * S(j) / (S(j)+E(j)+I(j)+R(j)) - b1 * lambda2(j) * (1-u1(j)) * S(j) / (S(j)+E(j)+I(j)+R(j)) + lambda * lambda2(j) + epsilon1 * lambda2(j) + mu * lambda2(j) - lambda * lambda3(j) - epsilon1 * lambda4(j) * (1+u2(j));
        n13 = -k2 - m3 * I(j) + b2 * lambda1(j) * (1-u1(j)) * S(j) / (S(j)+E(j)+I(j)+R(j)) - b2 * lambda2(j) * (1-u1(j)) * S(j) / (S(j)+E(j)+I(j)+R(j)) + epsilon2 * lambda3(j) * (1+u2(j)) + ( gamma + mu) - epsilon2 * lambda4(j);
        n14 = -m4 * R(j) + mu * lambda4(j);
        
        n21 = -m1 * 0.5*(S(j) + S(j-1)) + (lambda1(j)-h2*n11) * b1 * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5*(E(j) + E(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + (lambda1(j)-h2*n11) * b2 * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5*(I(j) + I(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + mu * (lambda1(j)-h2*n11) - b1 * (lambda2(j)-h2*n12) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 *(E(j) + E(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+ 0.5*(I(j)+I(j-1))+0.5*(R(j)+R(j-1))) - b2 * (lambda2(j)-h2*n12) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 *(I(j) + I(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+ 0.5*(I(j)+I(j-1))+0.5*(R(j)+R(j-1)));
        n22 = -k1 - m2 * 0.5*(E(j) + E(j-1)) + b1 * (lambda1(j)-h2*n11) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) - b1 * (lambda2(j)-h2*n12) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + lambda * (lambda2(j)-h2*n12) + epsilon1 * (lambda2(j)-h2*n12) + mu * (lambda2(j)-h2*n12) - lambda * (lambda3(j)-h2*n13) - epsilon1 * (lambda4(j)-h2*n14) * (1+(0.5*(u2(j) + u2(j-1))));
        n23 = -k2 - m3 * 0.5*(I(j) + I(j-1)) + b2 * (lambda1(j)-h2*n11) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) - b2 * (lambda2(j)-h2*n12) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + epsilon2 * (lambda3(j)-h2*n13) * (1+0.5*(u2(j)+u2(j-1))) + ( gamma + mu) - epsilon2 * (lambda4(j)-h2*n14);
        n24 = -m4 * 0.5*(R(j) + R(j-1)) + mu * (lambda4(j)-h2*n14);
        n31 = -m1 * 0.5*(S(j) + S(j-1)) + (lambda1(j)-h2*n21) * b1 * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5*(E(j) + E(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + (lambda1(j)-h2*n21) * b2 * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5*(I(j) + I(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + mu * (lambda1(j)-h2*n21) - b1 * (lambda2(j)-h2*n22) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 *(E(j) + E(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+ 0.5*(I(j)+I(j-1))+0.5*(R(j)+R(j-1))) - b2 * (lambda2(j)-h2*n22) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 *(I(j) + I(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+ 0.5*(I(j)+I(j-1))+0.5*(R(j)+R(j-1)));
        n32 = -k1 - m2 * 0.5*(E(j) + E(j-1)) + b1 * (lambda1(j)-h2*n21) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) - b1 * (lambda2(j)-h2*n22) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + lambda * (lambda2(j)-h2*n22) + epsilon1 * (lambda2(j)-h2*n22) + mu * (lambda2(j)-h2*n22) - lambda * (lambda3(j)-h2*n23) - epsilon1 * (lambda4(j)-h2*n24) * (1+(0.5*(u2(j) + u2(j-1))));
        n33 = -k2 - m3 * 0.5*(I(j) + I(j-1)) + b2 * (lambda1(j)-h2*n21) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) - b2 * (lambda2(j)-h2*n22) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + epsilon2 * (lambda3(j)-h2*n23) * (1+0.5*(u2(j)+u2(j-1))) + ( gamma + mu) - epsilon2 * (lambda4(j)-h2*n24);
        n34 = -m4 * 0.5*(R(j) + R(j-1)) + mu * (lambda4(j)-h2*n24);
        n41 = -m1 * 0.5*(S(j) + S(j-1)) + (lambda1(j)-h2*n31) * b1 * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5*(E(j) + E(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + (lambda1(j)-h2*n31) * b2 * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5*(I(j) + I(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + mu * (lambda1(j)-h2*n31) - b1 * (lambda2(j)-h2*n32) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 *(E(j) + E(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+ 0.5*(I(j)+I(j-1))+0.5*(R(j)+R(j-1))) - b2 * (lambda2(j)-h2*n32) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 *(I(j) + I(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+ 0.5*(I(j)+I(j-1))+0.5*(R(j)+R(j-1)));
        n42 = -k1 - m2 * 0.5*(E(j) + E(j-1)) + b1 * (lambda1(j)-h2*n31) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) - b1 * (lambda2(j)-h2*n32) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + lambda * (lambda2(j)-h2*n32) + epsilon1 * (lambda2(j)-h2*n32) + mu * (lambda2(j)-h2*n32) - lambda * (lambda3(j)-h2*n33) - epsilon1 * (lambda4(j)-h2*n34) * (1+(0.5*(u2(j) + u2(j-1))));
        n43 = -k2 - m3 * 0.5*(I(j) + I(j-1)) + b2 * (lambda1(j)-h2*n31) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) - b2 * (lambda2(j)-h2*n32) * (1-(0.5*(u1(j) + u1(j-1)))) * 0.5 * (S(j)+S(j-1)) / (0.5*(S(j)+S(j-1))+0.5*(E(j)+E(j-1))+0.5*(I(j)+I(j-1))+ 0.5*(R(j)+R(j-1))) + epsilon2 * (lambda3(j)-h2*n33) * (1+0.5*(u2(j)+u2(j-1))) + ( gamma + mu) - epsilon2 * (lambda4(j)-h2*n34);
        n44 = -m4 * 0.5*(R(j) + R(j-1)) + mu * (lambda4(j)-h2*n34);
        
        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);
    end
    
    % OPTIMALITY CONDITIONS
    U1 = max(0, min(1.0, ((lambda2 - lambda1).*b1.*S.*E + (lambda2 - lambda1).*b2.*S.*I ) ./ (l1 .* (S+E+I+R))));
    u1 = 0.01.*U1 + 0.99.*oldu1;
    U2 = max(0, min(1.0, ((lambda3 .* epsilon2 .* I - lambda4 .* epsilon1 .* E) ./(l2))));
    u2 = 0.01.*U2 + 0.99.*oldu2;
    
    % U1 = max(0, min(1, (L2 - L1).*b1.*x1.*x2./(b1)));
    % u1 = 0.01.*U1 + 0.99.*oldu1;
    % U2 = min(1.0, max(0, (((L2 - L3).*nu.*x2)./(b2))));
    % u2 = 0.01.*U2 + 0.99.*oldu2;
    % U3 = min(1.0, max(0, (((L1 - L7).*psi.*x1)./(b3))));
    % u3 = 0.01.*U3 + 0.99.*oldu3;
    
    % COST FUNCTION
    J     = k1*sum(E)*h + k2*sum(I)*h + l1*sum(u1.^2)*h + l2*sum(u2.^2)*h + m1./2*sum(S.^2)*h + m2./2*sum(E.^2)*h+ m3./2*sum(I.^2)*h+ m4/2*sum(R.^2)*h;
    Cost1 = k1.*cumsum(E)*h;                 % Total cost of threatened population
    Cost2 = k2.*cumsum(I)*h;
    Cost3 = l1.*cumsum(u1.^2)*h;
    Cost4 = l2.*cumsum(u2.^2)*h;
    Cost5 = m1./2.*cumsum(S.^2)*h;
    Cost6 = m2./2.*cumsum(E.^2)*h;                 % Total cost of control input u1
    Cost7 = m3./2.*cumsum(I.^2)*h;                 % Total cost of control input u2
    Cost8 = m4./2.*cumsum(R.^2)*h;                 % Total cost of control input u3
    J2    = Cost1 + Cost2 + Cost3 + Cost4 + Cost5 + Cost6 + Cost7+ Cost8;  % Cost at each time for ...plotting graphs
    %  % COST FUNCTION
    % J     = c1./2*sum(x4.^2)*h + b1./2*sum(u1.^2)*h + b2./2*sum(u2.^2)*h + b3./2*sum(u3.^2)*h + c2*max(x6);
    % Cost1 = c1./2.*cumsum(x4.^2)*h;                 % Total cost of threatened population
    % Cost2 = b1./2.*cumsum(u1.^2)*h;                 % Total cost of control input u1
    % Cost3 = b2./2.*cumsum(u2.^2)*h;                 % Total cost of control input u2
    % Cost4 = b2./2.*cumsum(u3.^2)*h;                 % Total cost of control input u3
    % Cost5 = c2.*x6;                                 % Total cost of deceased population
    % J2    = Cost1 + Cost2 + Cost3 + Cost4 + Cost5;  % Cost at each time for ...plotting graphs
    % 
    % CHECK CONVERGENCE TO STOP SWEEP METHOD
    temp1  = Delta*sum(abs(u1)) - sum(abs(oldu1 - u1));
    temp2  = Delta*sum(abs(u2)) - sum(abs(oldu2 - u2));
    temp3  = Delta*sum(abs(S)) - sum(abs(oldS - S));
    temp4  = Delta*sum(abs(E)) - sum(abs(oldE - E));
    temp5  = Delta*sum(abs(I)) - sum(abs(oldI - I));
    temp6  = Delta*sum(abs(R)) - sum(abs(oldR - R));
    temp7 = Delta*sum(abs(lambda1)) - sum(abs(oldlambda1 - lambda1));
    temp8 = Delta*sum(abs(lambda2)) - sum(abs(oldlambda2 - lambda2));
    temp9 = Delta*sum(abs(lambda3)) - sum(abs(oldlambda3 - lambda3));
    temp10 = Delta*sum(abs(lambda4)) - sum(abs(oldlambda4 - lambda4));
    test = min([temp1 temp2 temp3 temp4 temp5 temp6 temp7 temp8 temp9 temp10]);
end
disp(['number of loops: ' num2str(loopcnt)]);
disp(['Cost function: ' num2str(J)]);
disp(['Portion deceased: ' num2str(max(x6))]);
y(1,:)  = t;
y(2,:)  = S;
y(3,:)  = E;
y(4,:)  = I;
y(5,:)  = R;
y(9,:)  = lambda1;
y(10,:) = lambda2;
y(11,:) = lambda3;
y(12,:) = lambda4;
y(16,:) = u1;
y(17,:) = u2;
y(19,:) = J;
% IMMUNITY REACHED
% imm    = x5 + x7.*100;  % Percentage immune
% R_per = R*100;        % percentage threatened
% x6_per = x6*100;        % percentage deceased
plot(y(1,:), y(2,:)), grid on, xlabel('t'), ylabel('S')    % plotting S
plot(y(1,:), y(3,:)), grid on, xlabel('t'), ylabel('E')    % plotting E
plot(y(1,:), y(4,:)), grid on, xlabel('t'), ylabel('I')    % plotting I
plot(y(1,:), y(5,:)), grid on, xlabel('t'), ylabel('R')    % plotting R

I want the deterministic version of the code

[KALYAN ACHARJYA]

"I have a code for optimal control of a deterministic SEIR model. I am unable to code the stochastic version of the model" &

At end "I want the deterministic version of the code?"

Can you clarify more?

[priya anjali]

I have to implement noise terms in all the compartments S, E, I and R of the model. I have used Runge-Kutta method to solve the deterministic model(in the code provided by me). My question is where exactly should i add noise term in that code to get the plot for stochastic model?