Hello sir,
please check my code and resolve the error
mu = 0.012;
lambda = 1.2e-6;
omega = 0.75;
rho = 0.1;
K_T = 4e-4;
upsilon = 2.5e-4;
M_mt = 0.1;
alpha = 2.5e-3;
K_N = 4e-5;
epsilon = 1.2e-3;
beta = 0.3;
eta = 0.02;
sigma = 0.02;
K_D = 2e-4;
j = 4e-4;
gamma = 0.06;
P_LI = 2e-4;
K_L = 4e-4;
g = 1.2e-4;
v_l = 0.01;
% Time horizon
t0 = 0; tf = 12;
% Initial conditions
T0 = 1; N0 = 1e-5; D0 = 1e-5; L0 = 1e-5;
% Time delays
tau4 = 1; tau5 = 1; tau6 = 1; tau7 = 1;
% Solve the optimal control problem
x0 = [T0, N0, D0, L0];
[t, x] = ode45(@(t,x) ode_cancer(t, x, tau4, tau5, tau6, tau7, lambda, omega, rho, K_T, upsilon, M_mt, alpha, K_N, epsilon, beta, eta, sigma, K_D, j, gamma, P_LI, K_L, g, v_l), [t0, tf], x0);
% Extract the states
T = x(:,1);
N = x(:,2);
D = x(:,3);
L = x(:,4);
% Compute the optimal therapy function
Lambda = zeros(size(t));
for i = 1:length(t)
pT = -lambda*(1-mu*T(i)) + (omega*N(i) + rho*D(i) + gamma*L(i)) - K_T*Lambda(i)*T(i);
pN = upsilon + M_mt - (omega*T(i) - alpha*D(i))*N(i) - K_N*Lambda(i)*N(i) - epsilon*N(i);
pD = beta - (eta*L(i) + alpha*N(i) - sigma*T(i))*D(i) - K_D*Lambda(i)*D(i) - j*D(i);
pL = rho*D(i)*T(max(1,i-tau2))*T(max(1,i-tau4)) - gamma*L(i) - chi*N(i)*L(max(1,i-tau3))^2 + omega*N(i)*T(max(1,i-tau1))*T(max(1,i-tau4)) + P_LI - K_L*Lambda(i)*L(i) - g*L(i);
H = pT*Lambda(i) + pN*Lambda(i) + pD*Lambda(i) + pL*Lambda(i);
if H > 0
Lambda(i) = 1;
else
Lambda(i) = 0;
end
end
% Plot the results
figure;
subplot(2,2,1);
plot(t, T);
xlabel('Time (days)');
ylabel('Tumor Size (cm^3)');
title('Tumor Growth');
subplot(2,2,2);
plot(t, Lambda);
xlabel('Time (days)');
ylabel('Treatment');
title('Optimal Treatment');