graphs are not displaying

Hello so im trying to simulate a insulin infusion pump based on the stoljwik and hardy method, and i have this code running without an error, however both graphs are not displaying, what can i do to fix it.

% Parámetros del modelo de Stolwijk y Hardy
global k1 k2 k3 c theta alpha y betaa
k1 = 0.003;
k2 = 0.006;
k3 = 0.001;
c = 100;
theta = 2.5;
alpha = 7600;
y = 0;
betaa = 139000;
function main()
% Parámetros del controlador PID
Kp = 0.1;
Ki = 0.05;
Kd = 0.01;
int_error = 0;
prev_error = 0;
% Tiempo de simulación
tspan = [0 180];  % Desde 0 a 180 minutos
% Función de referencia (nivel deseado de glucosa)
reference = @(t) 100;
% Función de bomba de insulina
u = @(t, x) bomba_insulina(t, x);
% Simulación del sistema
[t, x] = ode45(@(t, x) modelo_stolwijk_hardy(t, x, u), tspan, initial_state);
% Resultados
glucosa = x(:, 1);
insulina = x(:, 2);
% Gráficos
clf();
figure;
subplot(211);
plot(t, glucosa);
xlabel('Tiempo (minutos)');
ylabel('Glucosa en sangre (mg/dL)');
title('Concentración de Glucosa en Sangre');
subplot(212);
plot(t, insulina);
xlabel('Tiempo (minutos)');
ylabel('Insulina en sangre (mU/L)');
title('Concentración de Insulina en Sangre');
% Mostrar las gráficas
grid on;
end
function dx = modelo_stolwijk_hardy(t, x, u)

G = x(1);
y = x(2);

% Utilizar la variable t
dx = [k1 * (u - G) - k2 * G;
    betaa * (G - theta) - alpha * y];
end
function u = bomba_insulina(t, x)
G = x(1);

% Condición de no producción de insulina
if G >= theta
    u = 0;
else
    u = betaa * (G - theta) - alpha * y;
end
end

Answer 1

[Torsten]

"initial_state" had to be defined.

globals had to be included where the variables were needed.

"y" has been removed from the globals because it seems to be a solution variable.

"u" has to be used as "u(t,x)", not simply as "u" in function "modelo_stolwijk_hardy".

"main" has to be called from the script part.

ode45 has been replaced by ode15s because your ODE system seems to be stiff.

"grid on" has been added to the first subplot.

I don't remember if maybe I forgot some more changes.

% Parámetros del modelo de Stolwijk y Hardy
global k1 k2 k3 c theta alpha betaa
k1 = 0.003;
k2 = 0.006;
k3 = 0.001;
c = 100;
theta = 2.5;
alpha = 7600;
betaa = 139000;
main()

function main()
% Parámetros del controlador PID
Kp = 0.1;
Ki = 0.05;
Kd = 0.01;
int_error = 0;
prev_error = 0;
% Tiempo de simulación
tspan = [0 180];  % Desde 0 a 180 minutos
% Función de referencia (nivel deseado de glucosa)
reference = @(t) 100;
% Función de bomba de insulina
u = @(t, x) bomba_insulina(t, x);
initial_state = [1;1];
% Simulación del sistema
[t, x] = ode15s(@(t, x) modelo_stolwijk_hardy(t, x, u), tspan, initial_state);
% Resultados
glucosa = x(:, 1);
insulina = x(:, 2);
% Gráficos
clf();
figure;
subplot(211);
plot(t, glucosa);
xlabel('Tiempo (minutos)');
ylabel('Glucosa en sangre (mg/dL)');
title('Concentración de Glucosa en Sangre');
grid on
subplot(212);
plot(t, insulina);
xlabel('Tiempo (minutos)');
ylabel('Insulina en sangre (mU/L)');
title('Concentración de Insulina en Sangre');
% Mostrar las gráficas
grid on;
end
function dx = modelo_stolwijk_hardy(t, x, u)
global k1 k2 k3 c theta alpha betaa
G = x(1);
y = x(2);

% Utilizar la variable t
dx = [k1 * (u(t,x) - G) - k2 * G;
    betaa * (G - theta) - alpha * y];
end
function u = bomba_insulina(t, x)
global k1 k2 k3 c theta alpha betaa
G = x(1);
y = x(2);

% Condición de no producción de insulina
if G >= theta
    u = 0;
else
    u = betaa * (G - theta) - alpha * y;
end
end