file_search_maxima : cons(sconcat("/home/renato/nuevos/divulgacion/coronavirus/numbers/modelo/###.{lisp,mac,mc}"), file_search_maxima)$

## Toy Model de una epidemia

kill(all)$

load(descriptive)$
lismax(li):=block([comp,mm],comp:length(li)-2,mm:1+sublist_indices(makelist(li[i] > li[i+1] or li[i+1] <= li[i+2],i,1,comp),lambda([x],x = 'false)))$
lismin(li):=block([comp,mm],comp:length(li)-2,mm:1+sublist_indices(makelist(li[i] < li[i+1] or li[i+1] >= li[i+2],i,1,comp),lambda([x],x = 'false)))$

### Modelo básico 

kill(alpha,beta,delta,gamma,s,i,r,f)$

Ecuaciones del modelo básico

eq1:-alpha*s*i/(s+i+r)+delta*r;
eq2:alpha*s*i/(s+i+r)-(beta)*i;
eq3:beta*i-delta*r;

eq1+eq2+eq3,ratsimp;

Soluciones estacionarias

esta:solve([eq1=0,eq2=0,eq3=0],[s,i,r]);

/* Solucion estacionaria */
esta:solve([eq1=0,eq2=0,eq3=0],[s,i,r])$
expand(esta[1]);
ests1:rhs(esta[1][1]);esti1:rhs(esta[1][2]);estr1:rhs(esta[1][3]);
expand(esta[2]);
ests2:rhs(esta[2][1]);esti2:rhs(esta[2][2]);estr2:rhs(esta[2][3]);

Análisis de estabilidad de las soluciones estacionarias

define(se2(alpha,beta,delta,N0), (N0*beta)/(alpha) );

define(ie2(alpha,beta,delta,N0), -(N0*(beta-alpha)*delta)/((delta+beta)*alpha) );

define(re2(alpha,beta,delta,N0), -(N0*(beta-alpha)*beta)/((delta+beta)*alpha));

Matriz jacobiana del sistema

jac:jacobian([eq1,eq2,eq3],[s,i,r]);

Estabilidad de la primera solución estacionaria

mat1:ev(jac,[s=ests1,i=esti1,r=estr1] );
emat1:eigenvalues(%);

ests1+esti1+estr1;

define(autovalores1(alpha,beta,delta), emat1[1]);

Estabilidad de la segunda solución estacionaria

mat2:ev(jac,[s=ests2,i=esti2,r=estr2] );
emat2:eigenvalues(%);

ests2+esti2+estr2,ratsimp;

Autovalores de la segunda solución estacionaria en función de las variables alpha,beta,delta

define(autovalores2(alpha,beta,delta), emat2[1]);

delta^4+(4·beta-2·alpha)·delta^3+(8·beta^2-8·alpha·beta+alpha^2)·delta^2+(4·beta^3-4·alpha·beta^2)·delta;

Del:expand(%);

Del-delta^2*(alpha+delta)^2,expand;

factor(%);

## Simulación

define(fs(t,alpha,delta,s,i,r),eq1);  /* susceptibles */
define(fi(t,alpha,beta,gamma,s,i),eq2);  /* infectados */
define(fr(t,beta,delta,i,r), eq3 ); /* recuperados */ 
define(ff(t,gamma,i),gamma*i) ; /* fallecidos */

fs(t,alpha,delta,s,i,r);fi(t,alpha,beta,gamma,s,i);fr(t,beta,delta,i,r);

### Epidemia éndemica

Valores de las constantes del modelo

alpha:0.065; beta:0.03; delta:0.02; gamma:0.02; 

Si alpha > beta nos encontramos con la segunda solución estacionaria y en el caso constrario con la primera

alpha-beta;

realpart(autovalores1(alpha,beta,delta));

realpart(autovalores2(alpha,beta,delta));

Valores iniciales de los diferentes bloques

s0:10000$ i0:100$ r0:0$ f0:0$ /* Initial values */;
N0:s0+i0+f0;

steo:se2(alpha,beta,delta,N0);
iteo:ie2(alpha,beta,delta,N0);
rteo:re2(alpha,beta,delta,N0);

se2(alpha,beta,delta,N0)+ie2(alpha,beta,delta,N0)+re2(alpha,beta,delta,N0);

Resolución del sistema diferencial

sol:rk( [fs(t,alpha,delta,s,i,r),fi(t,alpha,beta,gamma,s,i),fr(t,beta,delta,i,r)] ,  
    [s,i,r],[s0,i0,r0],[t,0,700,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0);

tt0[1401];

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data1.dat" )$

concat("susceptibles inicio=",s0,"; fin=", ss:float( ceiling(1000*ss0[lll])/1000 ), "; teo=",steo  );
concat("enfermas inicio=",i0,"; fin=",ii:float( ceiling(1000*ii0[lll])/1000 ), "; teo=",iteo);
concat("recuperados inicio=",r0,"; fin=", rr:float( ceiling(100*rr0[lll])/100 ), "; teo=",rteo);
/*  concat("fallecidos inicio=",f0,"; fin=", ff:float( ceiling(1000*ff0[lll])/1000 ), "; teo=",fteo); */
float( floor( 100*(ss+ii+rr)) /100);

Mf(lista,p0,pf):= block([ll], ll:makelist(lista[k],k,p0,pf), apply ("+", ll)/length(ll) ) ;

ini:301;fin: 1401;

tt0[ini]; tt0[fin]; 

Mf(ff0,ini,fin);

escala:max(smax(ss0))$
g1a:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", yrange = [0, escala*1.2],  
            /* title=concat("susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , */
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0),  line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$

wxdraw(g1a)$

draw(g1a)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "basic-asint")$

escala:max(smax(ss0))$
escalai:smax(ii0)$
escalaf:smax(ii0)*gamma$

draw2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",   yrange = [0, escalaf*1.2],  
                      point_type = filled_circle, point_size =0.5, 
                       color =  magenta,   key="fallecidos",  points(tt0,ff0) )$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "basic-asintf")$

draw2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",   yrange = [0, escalaf*1.2],  
                      point_type = filled_circle, point_size =0.5, 
                       color =  magenta,   key="fallecidos",  points(tt0,ff0) )$

### Epidemia extinguida

Valores de las constantes

alpha:0.065; beta:0.07; delta:0.02; gamma:0.02; 

alpha-beta;

realpart(autovalores1(alpha,beta,delta));

realpart(autovalores2(alpha,beta,delta));

s0:10000$ i0:500$ r0:0$ f0:0$ /* Initial values */;
N0:s0+i0+f0;

steo:N0;
iteo:0;
rteo:0;

sol:rk( [fs(t,alpha,delta,s,i,r),fi(t,alpha,beta,gamma,s,i),fr(t,beta,delta,i,r)] ,  
    [s,i,r],[s0,i0,r0],[t,0,700,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][4],k,1,length(sol),5)$ 
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data1a.dat" )$

concat("susceptibles inicio=",s0,"; fin=", ss:float( ceiling(1000*ss0[lll])/1000 ), "; teo=",steo  );
concat("enfermas inicio=",i0,"; fin=",ii:float( ceiling(1000*ii0[lll])/1000 ), "; teo=",iteo);
concat("recuperados inicio=",r0,"; fin=", rr:float( ceiling(100*rr0[lll])/100 ), "; teo=",rteo);
/*  concat("fallecidos inicio=",f0,"; fin=", ff:float( ceiling(1000*ff0[lll])/1000 ), "; teo=",fteo); */
float( floor( 100*(ss+ii+rr)) /100);

s0+i0+r0;
ss+ii+rr;

escala:max(smax(ss0))$
g2:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",  yrange = [0, escala*1.35],  
            title=concat("susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(g2)$

draw(g2)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "basicno")$

## Confinamiento

kill(all)$
load(descriptive)$

Función confinamiento

phi(t,alpha_1, alpha_2, ti, tf,tau) := if  (t>=ti and t<ti + tau) then (alpha_2-alpha_1)*(t-ti)/(tau) + alpha_1
	else if (t>=ti + tau and t<tf) then alpha_2 else if (t>=tf and t<tf+tau) then (alpha_1-alpha_2)*(t-tf)/(tau) + alpha_2 
	else alpha_1;
Mf(lista,p0,pf):= block([ll], ll:makelist(lista[k],k,p0,pf), apply ("+", ll)/length(ll) ) ;

fs(t,alpha, delta,s,i,r) := -alpha*s*i/(s+i+r)+delta*r;  /* susceptibles */
fi(t,alpha,beta,gamma,s,i) := alpha*s*i/(s+i+r)-(beta)*i;  /* sanas */
fr(t,beta,delta,i,r) := beta*i-delta*r;  /* enfermas */ 
ff(t,gamma,i):= gamma*i; /* inmunes */

### Confinamiento permanente: sin acabar la epidemia

fs(t,alpha, delta,s,i,r) ;
fi(t,alpha,beta,gamma,s,i) ;
fr(t,beta,delta,i,r) := beta*i-delta*r;    

alpha_1:0.065$ alpha_2 : 0.05$ beta:0.03$ delta:0.02$ gamma:0.02$  tau:50$ ti:150$  tf:700$ /* confinamiento permanente*/ 
s0:10000$ i0:100$ r0:0$ f0:0$ /* Initial values */;

phi(t, alpha_1, alpha_2, ti, tf,tau);

sol:rk( [fs(t,phi(t, alpha_1, alpha_2, ti, tf,tau),delta,s,i,r),fi(t,phi(t, alpha_1, alpha_2, ti, tf,tau),beta,gamma,s,i),fr(t,beta,delta,i,r)] ,  
    [s,i,r],[s0,i0,r0],[t,0,700,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

iper:ii0$
fper:ff0$

tt0[1401];

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data2.dat" )$

Gráfica confinamiento permanente

escala:max(smax(ss0))$
gcp:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",  yrange = [0, escala*1.35],  
            title=concat("susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(gcp)$

ini:301;fin: 1401;

tt0[ini]; tt0[fin]; 

Mf(iper,ini,fin);

Mf(fper,ini,fin);

draw(gcp)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "confino")$

### Confinamiento permanente: eliminando la epidemia

alpha_1:0.065$ alpha_2 : 0.025$ beta:0.03$ delta:0.02$ gamma:0.02$  tau:50$ ti:150$  tf:700$ /* confinamiento permanente*/ 
s0:10000$ i0:100$ r0:0$ f0:0$ /* Initial values */;

sol:rk( [fs(t,phi(t, alpha_1, alpha_2, ti, tf,tau),delta,s,i,r),fi(t,phi(t, alpha_1, alpha_2, ti, tf,tau),beta,gamma,s,i),fr(t,beta,delta,i,r)] ,  
    [s,i,r],[s0,i0,r0],[t,0,700,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data2a.dat" )$

Gráfica confinamiento permanente, con desaparición del virus

escala:max(smax(ss0))$
gcp:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",  yrange = [0, escala*1.35],  
            title=concat("susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(gcp)$

draw(gcp)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "confisi")$

phi(t, alpha_1, alpha_2, ti, tf,tau);

Gráfica de la función confinamiento

wxdraw2d(line_width  = 3,  explicit(phi(t, 1, 1/2, 50,200,50),t,0,300),  yrange = [.4, 1.1]   );

draw2d(line_width  = 3,  explicit(phi(t, 1, 1/2, 50,200,50),t,0,300),  yrange = [.4, 1.1]   )$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "alpha_t")$

### Cofinamiento parcial: olas de contagios

alpha_1:0.065$ alpha_2 : 0.05$ beta:0.03$ delta:0.02$ gamma:0.02$ tau : 50$ ti : 150$  tf:300$ /* confinamiento parcial*/
s0:10000$ i0:100$ r0:0$ f0:0$ /* Initial values */;

sol:rk( [fs(t,phi(t, alpha_1, alpha_2, ti, tf,tau),delta,s,i,r),fi(t,phi(t, alpha_1, alpha_2, ti, tf,tau),beta,gamma,s,i),fr(t,beta,delta,i,r)] ,  
    [s,i,r],[s0,i0,r0],[t,0,700,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

i150:ii0$
f150:ff0$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data3.dat" )$

Gráfica del confinamiento parcial con finalización del confinamiento en 300

escala:max(smax(ss0))$
gcp:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",  yrange = [0, escala*1.35],  
            title=concat("susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(gcp)$

draw(gcp)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "confi150")$

alpha_1:0.065$ alpha_2 : 0.05$beta:0.03$ delta:0.02$ gamma:0.02$ tau : 50$ ti : 150$  tf:450$ /* confinamiento parcial*/
s0:10000$ i0:100$ r0:0$ f0:0$ /* Initial values */;

sol:rk( [fs(t,phi(t, alpha_1, alpha_2, ti, tf,tau),delta,s,i,r),fi(t,phi(t, alpha_1, alpha_2, ti, tf,tau),beta,gamma,s,i),fr(t,beta,delta,i,r)] ,  
    [s,i,r],[s0,i0,r0],[t,0,700,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

i300:ii0$
f300:ff0$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data3a.dat" )$

Gráfica del confinamiento parcial con finalización del confinamiento en 450

escala:max(smax(ss0))$
gcp:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",  yrange = [0, escala*1.35],  
            title=concat("susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(gcp)$

draw(gcp)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "confi300")$

Comparación de los infectados en los casos anteriores

gcompi:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",  yrange = [0, escala*0.3],  
            point_type = filled_circle, point_size =.5, 
            color = blue, key="infectados 150", points(tt0,i150), 
             color = red,   key="infectados 300",  points(tt0,i300) )$
wxdraw(gcompi);

Comparación de los fallecidos en los casos anteriores

gcompf:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",  yrange = [0, 60],  
            point_type = filled_circle, point_size =.5, 
            color = blue, key="fallecidos 150", points(tt0,f150), 
             color = red,   key="fallecidos 300",  points(tt0,f300),
             color = black,   key="fallecidos ∞",  points(tt0,fper) )$
wxdraw(gcompf);

length(fper);

draw(gcompf)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "conficomp")$

Mf(lista,p0,pf):= block([ll], ll:makelist(lista[k],k,p0,pf), apply ("+", ll)/length(ll) ) ;

tt0[2000/5+1]; tt0[7000/5+1];

Mf(f150,ini,fin);
Mf(f300,ini,fin);
Mf(fper,ini,fin);

Mf(i150,ini,fin);
Mf(i300,ini,fin);
Mf(iper,ini,fin);

Mf(fper,ini,fin)/Mf(f150,ini,fin);

Mf(iper,ini,fin)/Mf(i150,ini,fin);

### Confinamiento estricto no indefinido

alpha_1:0.065$ alpha_2 : 0.025$beta:0.03$ delta:0.02$ gamma:0.02$ tau : 50$ ti : 150$  tf:450$ /* confinamiento parcial*/
s0:10000$ i0:100$ r0:0$ f0:0$ /* Initial values */;

sol:rk( [fs(t,phi(t, alpha_1, alpha_2, ti, tf,tau),delta,s,i,r),fi(t,phi(t, alpha_1, alpha_2, ti, tf,tau),beta,gamma,s,i),fr(t,beta,delta,i,r)] ,  
    [s,i,r],[s0,i0,r0],[t,0,700,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

in300:ii0$
fn300:ff0$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data4.dat" )$

Confinamiento parcial estricto

escala:max(smax(ss0))$
gcp:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",  yrange = [0, escala*1.35],  
            title=concat("susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(gcp)$

draw(gcp)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "conficasi")$

kill(leve,fuerte);

ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50)) $
leve:append([lev] ,makelist(i300[k],k,1,length(i300),10)) $
fuerte:append([fue],makelist(in300[k],k,1,length(in300),10)) $
write_data (makelist([ttl[k],leve[k],fuerte[k]],k,1,length(ttl)), "data4a.dat" )$

Comparación infectados con medidas leves y con medidas estrictas

gcompi:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",  yrange = [0, escala*0.3],  
            point_type = filled_circle, point_size =.5, 
            color = blue, key="infectados medidas leves", points(tt0,i300), 
             color = red,   key="infectados medidas fuertes",  points(tt0,in300) )$
wxdraw(gcompi);

write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data3.dat" )$

draw(gcompi);
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "confifuerte")$

## Función Infección

load("descriptive")$

Ecuaciones del sistema

fs(t,alpha, delta,s,i,r) := -alpha*s*i/(s+i+r)+delta*r;  /* susceptibles */
fi(t,alpha,beta,s,i) := alpha*s*i/(s+i+r)-(beta)*i;  /* sanas */
fr(t,beta,delta,i,r) := beta*i-delta*r;  /* enfermas */ 
ff(t,gamma,i):= gamma*i; /* inmunes */

Definición de la función Confinamiento y la función mortalidad

phi(t,alpha_1, alpha_2, ti, tf,tau) := if  (t>=ti and t<ti + tau) then (alpha_2-alpha_1)*(t-ti)/(tau) + alpha_1
	else if (t>=ti + tau and t<tf) then alpha_2 else if (t>=tf and t<tf+tau) then (alpha_1-alpha_2)*(t-tf)/(tau) + alpha_2 
	else alpha_1;
Mf(lista):= block([ll], ll:makelist(lista[k],k,1,2000), apply ("+", ll)/length(ll) ) ;

Definición de la función infección

psi(alpha,  kappa, phip) :=alpha*(1-c*phip)^kappa;

Definición de las funciones φ_1 y φ_2

phip1(x) :=if x <= 0.1 then 0 else if (0.1 < x and x <= 0.3) then 1/5 else if (0.3 < x and x <= 0.5) then 2/5 else 3/5; /*φ_1 escalonada*/

phip2(x) := x; /*φ_2 lineal*/

wxdraw2d(grid = true,  line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
        xlabel="% Infectados",   title=concat("Funcion φ_1 y φ_2") , 
          explicit(phip1(x), x, 0,1), color=red,
    explicit(phip2(x), x, 0,1), yrange=[-0.1, 1]);

draw2d(grid = true,  line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
        xlabel="% Infectados",   title=concat("Funcion φ_1 y φ_2)") , yrange=[-0.1, 1],
          explicit(phip1(x), x, 0,1), color=red,
          explicit(phip2(x), x, 0,1));
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "phip")$

### Ejemplos

### φ, κ=0 (comportamiento negligente)

alpha_1:0.065$ beta:0.03$ delta:0.02$ gamma:0.02$ c: 0.01$ kappa:0$
s0:10000$ i0:100$ r0:0$ f0:0$  /* Initial values */;

sol:rk( [fs(t,psi(alpha_1, kappa, phip1(i/(s+i+r))),delta,s,i,r),
            fi(t,psi(alpha_1, kappa, phip1(i/(s+i+r))),beta,s,i),
            fr(t,beta,delta,i,r)],   [s,i,r],[s0,i0,r0],[t,0,500,0.05] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

Mf(ii0*gamma);

i00:makelist(sol[k][3],k,1,length(sol),50)$
i00n:append([inf00] ,%)$

escala:max(smax(ss0))$
gk0:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( "κ=", kappa,  ", susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(gk0);

draw(gk0)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "gk0")$

### φ escalonada, κ=100 (comportamiento moderado)

alpha_1:0.065$ beta:0.03$ delta:0.02$ gamma:0.02$ c: 0.01$ kappa:100$
s0:10000$ i0:100$ r0:0$ f0:0$  /* Initial values */;

sol:rk( [fs(t,psi(alpha_1, kappa, phip1(i/(s+i+r))),delta,s,i,r),
            fi(t,psi(alpha_1, kappa, phip1(i/(s+i+r))),beta,s,i),
            fr(t,beta,delta,i,r)],   [s,i,r],[s0,i0,r0],[t,0,500,0.05] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data5.dat" )$

i100:makelist(sol[k][3],k,1,length(sol),50)$
i100n:append([inf100] ,%)$

Mf(ii0*gamma);

escala:max(smax(ss0))$
gek100:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( "κ=", kappa,  ", susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(gek100);

draw(gek100)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "gek100")$

### φ escalonada, κ=200 (comportamiento bueno)

alpha_1:0.065$ beta:0.03$ delta:0.02$ gamma:0.02$ c: 0.01$ kappa:200$
s0:10000$ i0:100$ r0:0$ f0:0$  /* Initial values */;

sol:rk( [fs(t,psi(alpha_1, kappa, phip1(i/(s+i+r))),delta,s,i,r),
            fi(t,psi(alpha_1, kappa, phip1(i/(s+i+r))),beta,s,i),
            fr(t,beta,delta,i,r)],   [s,i,r],[s0,i0,r0],[t,0,500,0.05] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data5a.dat" )$

i200:makelist(sol[k][3],k,1,length(sol),50)$
i200n:append([inf200] ,%)$

Mf(ff0);

escala:max(smax(ss0))$
gek200:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( "κ=", kappa,  ", susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(gek200);

draw(gek200)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "gek200")$

### φ lineal, κ=100 (comportamiento moderado)

alpha_1:0.065$ beta:0.03$ delta:0.02$ gamma:0.02$ c: 0.01$ kappa:100$
s0:10000$ i0:100$ r0:0$ f0:0$  /* Initial values */;

sol:rk( [fs(t,psi(alpha_1, kappa, phip2(i/(s+i+r))),delta,s,i,r),
            fi(t,psi(alpha_1, kappa, phip2(i/(s+i+r))),beta,s,i),
            fr(t,beta,delta,i,r)],   [s,i,r],[s0,i0,r0],[t,0,500,0.05] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data5b.dat" )$

il100:iil$

il100:makelist(sol[k][3],k,1,length(sol),50)$
il100n:append([infil100] ,%)$

Mf(ii0*gamma);

escala:max(smax(ss0))$
glk100:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( "κ=", kappa,  ", susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(glk100);

draw(glk100)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "glk100")$

### φ lineal, κ=200 (comportamiento bueno)

alpha_1:0.065$ beta:0.03$ delta:0.02$ gamma:0.02$ c: 0.01$ kappa:200$
s0:10000$ i0:100$ r0:0$ f0:0$  /* Initial values */;

sol:rk( [fs(t,psi(alpha_1, kappa, phip2(i/(s+i+r))),delta,s,i,r),
            fi(t,psi(alpha_1, kappa, phip2(i/(s+i+r))),beta,s,i),
            fr(t,beta,delta,i,r)],   [s,i,r],[s0,i0,r0],[t,0,500,0.05] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data5c.dat" )$

il200:makelist(sol[k][3],k,1,length(sol),50)$
il200n:append([infil200] ,%)$

ttn:makelist(sol[k][1],k,1,length(sol),50)$

Mf(ii0)*gamma;

escala:max(smax(ss0))$
glk200:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( "κ=", kappa,  ", susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(glk200);

draw(glk200)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "glk200")$

ratprint: false$

escala:max(smax(i100),smax(i200),smax(i00))$
compa3:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
        xlabel="unidades de tiempo",  
        title=concat( "Evolución de los enfermos para distintos valores de k") , 
        point_type = filled_circle, point_size =0.5, 
        color = red, key="k=0", points(ttn,i00),        
        color = blue, key="lineal, k=100", points(ttn,il100),  
        color = black, key="escalonado, k=100", points(ttn,i100) )$
wxdraw(%);

length(i00);

escala:max(smax(i100),smax(i200),smax(i00))$
compa3:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",   
            title=concat( "Evolución de los enfermos para distintos valores de k") , 
            point_type = filled_circle, point_size =0.5, 
            color = red, key="k=0", points(ttn,i00),
            color = blue, key="lineal, k=200", points(ttn,il200),   
            color = black, key="escalonado, k=200", points(ttn,i200) )$
wxdraw(%);

### φ lineal, κ = I/10 (κ dinámica)

alpha_1:0.065$ beta:0.03$ delta:0.02$ gamma:0.02$ c: 0.01$ kappa:i/10$
s0:10000$ i0:100$ r0:0$ f0:0$  /* Initial values */;

sol:rk( [fs(t,psi(alpha_1, kappa, phip2(i/(s+i+r))),delta,s,i,r),
            fi(t,psi(alpha_1, kappa, phip2(i/(s+i+r))),beta,s,i),
            fr(t,beta,delta,i,r)],   [s,i,r],[s0,i0,r0],[t,0,500,0.05] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data5d.dat" )$

ilki:makelist(sol[k][3],k,1,length(sol),50)$
ilkin:append([infilki] ,%)$

Mf(ii0*gamma);

escala:max(smax(ss0))$
glki:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( "κ=I/10",  ", susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(glki);

draw(glki)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "glki")$

write_data (makelist([ttl[k],i00n[k],il100n[k],ilkin[k],il200n[k] ],k,1,length(ttl)), "data5e.dat" )$

escala:max(smax(ilki),smax(i00))$
compa4:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",   
            title=concat( "Evolución de los enfermos para distintos valores de k") , 
            point_type = filled_circle, point_size =0.5, 
            color = red, key="k=0", points(ttn,i00), 
            color = blue, key="k=100", points(ttn,il100) ,
            color = dark_green , key="k=I/10", points(ttn,ilki),  
            color =black, key="k=200", points(ttn,il200) 
)$
wxdraw(%);

draw(compa4)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "compa4")$

### Caso con confinamiento permanente φ lineal, κ = 200 (comportamiento responsable)

alpha_1:0.065$ alpha_2 : 0.05$ beta:0.03$ delta:0.02$ gamma:0.02$  tau:50$ ti:150$  tf:2000$ c: 0.01$ kappa:200 $
s0:10000$ i0:100$ r0:0$ f0:0$ /* Initial values */;

sol:rk( [fs(t,psi(phi(t, alpha_1, alpha_2, ti, tf,tau), kappa, phip2(i/(s+i+r))),delta,s,i,r),
          fi(t,psi(phi(t, alpha_1, alpha_2, ti, tf,tau), kappa, phip2(i/(s+i+r))),beta,s,i),
          fr(t,beta,delta,i,r)] ,      [s,i,r],[s0,i0,r0],[t,0,500,0.05] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data5f.dat" )$

ipk:makelist(sol[k][3],k,1,length(sol),50)$
ipkn:append([vipk] ,%)$

Mf(ii0)*gamma;

escala:max(smax(ss0))$
gcpk200:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",    
            title=concat("susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(gcpk200);

draw(gcpk200)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "gcpk200")$

### Caso con confinamiento permanente φ, κ = 0 (comportamiento negligente)

alpha_1:0.065$ alpha_2 : 0.05$ beta:0.03$ delta:0.02$ gamma:0.02$  tau:50$ ti:150$  tf:2000$ c: 0.01$ kappa:0 $
s0:10000$ i0:100$ r0:0$ f0:0$ /* Initial values */;

sol:rk( [fs(t,psi(phi(t, alpha_1, alpha_2, ti, tf,tau), kappa, phip2(i/(s+i+r))),delta,s,i,r),
          fi(t,psi(phi(t, alpha_1, alpha_2, ti, tf,tau), kappa, phip2(i/(s+i+r))),beta,s,i),
          fr(t,beta,delta,i,r)] ,  
    [s,i,r],[s0,i0,r0],[t,0,500,0.05] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "data5d.dat" )$

ipk0:makelist(sol[k][3],k,1,length(sol),50)$
ipk0n:append([vipk0] ,%)$

Mf(ii0)*gamma;

escala:max(smax(ss0))$
gcpk0:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",    
            title=concat("susceptibles (azul),", " enfermos  (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(gcpk0);

draw(gcpk0)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "gcpk0")$

write_data (makelist([ttl[k],ipkn[k],ipk0n[k]],k,1,length(ttl)), "data5g.dat" )$

escala:max(smax(ipk0),smax(ipk))$
compa5:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo",    
            title=concat( "Evolución de los enfermos para distintos valores de k con confinamiento") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="κ=200", points(ttn,ipk),   
                      color = red, key="κ=0", points(ttn,ipk0) )$
wxdraw(%);

draw(compa5)$
draw_file(dimensions = [1000,700],  terminal = 'pdf, file_name = "compa5")$

## Vacunas

load("descriptive")$
ratprint: false$

Ecuaciones del sistema

fs(t,alpha, delta,s,i,r) := -alpha*s*i/(s+i+r)+delta*r;  /* susceptibles */
fi(t,alpha,beta,s,i) := alpha*s*i/(s+i+r)-(beta)*i;  /* sanas */
fr(t,beta,delta,i,r) := beta*i-delta*r;  /* enfermas */ 
ff(t,gamma,i):= gamma*i; /* inmunes */

Mf(lista,p0,pf):= block([ll], ll:makelist(lista[k],k,p0,pf), apply ("+", ll)/length(ll) ) ;

ini:301;fin: 1401;

Existen diferentes vías para la implementación de la vacunación en el modelo

- disminución de la tasa de contagio
- aumento de la tasa de recuperación
- eliminación de personas susceptibles

### Disminución de la tasa de contagio y aumento de la de recuperación

En este caso crearemos una nueva función que es lineal y que depende de la velocidad (tiempo en el que toda la población estaria vacunada) de vacunación eta, es similar a la funcion confinamiento pero aplicada a un valor mayor en el tiempo

η_α(t, tv, nu, alpha_1, alpha_2) := if t< tv then alpha_1 else max((alpha_2-alpha_1)/(nu)*(t-tv)+alpha_1, alpha_2);

alpha_1:0.065$ alpha_2 : 0.05$beta_1:0.03$  beta_2 : 0.04$  delta:0.02$ gamma:0.02$ c: 0.01$ kappa:0$ nu : 400 $ tv:200 $
s0:10000$ i0:100$ r0:0$ f0:0$  /* Initial values */;

wxdraw2d(grid = true,  line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
        xlabel="unidades de tiempo",   title=concat("Funcion α(t)") , 
         explicit(η_α(t, tv, 400, alpha_1, alpha_2), t, 0, 800),  
    color=red,explicit(η_α(t, tv, 100, alpha_1, alpha_2), t, 0, 800),
   yrange=[alpha_2-0.005, alpha_1+0.005]);

En este caso crearemos una nueva función que es lineal y que depende de la velocidad (tiempo en el que toda la población estaria vacunada) de vacunación nu

σ_β(t, tv, nu, beta_1, beta_2) := if t< tv then beta_1 else min((beta_2-beta_1)/(nu)*(t-tv)+beta_1, beta_2);

wxdraw2d(grid = true,  line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
        xlabel="unidades de tiempo",   title=concat("Funcion β(t)") , 
         explicit(σ_β(t, tv, 400, beta_1, beta_2), t, 0, 800),
color=red,    explicit(σ_β(t, tv, 100, beta_1, beta_2), t, 0, 800),
   yrange=[beta_1-0.005, beta_2+0.005]);

### Modelo básico (sin vacunas)

alpha_1:0.065$ alpha_2 : 0.05$beta_1:0.03$  beta_2 : 0.04$  delta:0.02$ gamma:0.02$ c: 0.01$ kappa:0$ nu : 400 $ tv:200 $
s0:10000$ i0:100$ r0:0$ f0:0$  /* Initial values */;

sol:rk( [fs(t,alpha_1,delta,s,i,r),
            fi(t,alpha_1,beta_1,s,i),
            fr(t,beta_1,delta,i,r)],   [s,i,r],[s0,i0,r0],[t,0,800,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),50)$ 
ii0:makelist(sol[k][3],k,1,length(sol),50)$ 
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */

ttl0:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
iil0:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$

length(makelist(sol[k][3],k,1,length(sol),5));

tt0[1401];

Mf(makelist(sol[k][3],k,1,length(sol),5),ini,fin)*gamma;

### Modelo híbrido α(t), β(t)

### Vacuna efectiva

alpha_1:0.065$alpha_2 : 0.05$ beta_1:0.03$ beta_2 : 0.04 $ delta:0.02$ gamma:0.02$ c: 0.01$ kappa:0$nu : 400 $ tv:200 $
s0:10000$ i0:100$ r0:0$ f0:0$  /* Initial values */;

σ_β(t, tv, nu, beta_1, beta_2) := if t< tv then beta_1 else min((beta_2-beta_1)/(nu)*(t-tv)+beta_1, beta_2);
η_α(t, tv, nu, alpha_1, alpha_2) := if t< tv then alpha_1 else max((alpha_2-alpha_1)/(nu)*(t-tv)+alpha_1, alpha_2);

sol:rk( [fs(t,η_α(t, tv, nu, alpha_1, alpha_2),delta,s,i,r),
            fi(t,η_α(t, tv, nu, alpha_1, alpha_2),σ_β(t, tv, nu, beta_1, beta_2),s,i),
            fr(t,σ_β(t, tv, nu, beta_1, beta_2),delta,i,r)],   [s,i,r],[s0,i0,r0],[t,0,800,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "datav1.dat" )$

ve:makelist(sol[k][3],k,1,length(sol),50)$
ven:append([vace] ,%)$

gk04:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( " susceptibles (azul),", " enfermos (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(gk04);

Mf(ii0,ini,fin)*gamma;

### Vacuna ineficaz 

alpha_1:0.065$alpha_2 : 0.06$ beta_1:0.03$ beta_2 : 0.035 $ delta:0.02$ gamma:0.02$ c: 0.01$ kappa:0$nu : 400 $ tv:200 $
s0:10000$ i0:100$ r0:0$ f0:0$  /* Initial values */;

σ_β(t, tv, nu, beta_1, beta_2) := if t< tv then beta_1 else min((beta_2-beta_1)/(nu)*(t-tv)+beta_1, beta_2);
η_α(t, tv, nu, alpha_1, alpha_2) := if t< tv then alpha_1 else max((alpha_2-alpha_1)/(nu)*(t-tv)+alpha_1, alpha_2);

sol:rk( [fs(t,η_α(t, tv, nu, alpha_1, alpha_2),delta,s,i,r),
            fi(t,η_α(t, tv, nu, alpha_1, alpha_2),σ_β(t, tv, nu, beta_1, beta_2),s,i),
            fr(t,σ_β(t, tv, nu, beta_1, beta_2),delta,i,r)],   [s,i,r],[s0,i0,r0],[t,0,800,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

vi:makelist(sol[k][3],k,1,length(sol),50)$
vin:append([vaci] ,%)$

gk06:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( " susceptibles (azul),", " enfermos (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(gk06);

Mf(ii0,ini,fin)*gamma;

### Vacuna totalmente efectiva - fin epidemia

alpha_1:0.065$alpha_2 : 0.04$ beta_1:0.03$ beta_2 : 0.05 $ delta:0.02$ gamma:0.02$ c: 0.01$ kappa:0$nu : 400 $ tv:200 $
s0:10000$ i0:100$ r0:0$ f0:0$  /* Initial values */;

σ_β(t, tv, nu, beta_1, beta_2) := if t< tv then beta_1 else min((beta_2-beta_1)/(nu)*(t-tv)+beta_1, beta_2);
η_α(t, tv, nu, alpha_1, alpha_2) := if t< tv then alpha_1 else max((alpha_2-alpha_1)/(nu)*(t-tv)+alpha_1, alpha_2);

sol:rk( [fs(t,η_α(t, tv, nu, alpha_1, alpha_2),delta,s,i,r),
            fi(t,η_α(t, tv, nu, alpha_1, alpha_2),σ_β(t, tv, nu, beta_1, beta_2),s,i),
            fr(t,σ_β(t, tv, nu, beta_1, beta_2),delta,i,r)],   [s,i,r],[s0,i0,r0],[t,0,800,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$ 
rr0:makelist(sol[k][4],k,1,length(sol),5)$ 
ff0:makelist(gamma*sol[k][3],k,1,length(sol),5)$ 
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
ffl:append([fal] , makelist(gamma*sol[k][3],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],ffl[k] ],k,1,length(ttl)), "datav2.dat" )$

vme:makelist(sol[k][3],k,1,length(sol),50)$
vmen:append([vacme] ,%)$

gk07:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( " susceptibles (azul),", " enfermos (rojo),",  " recuperados (negro)") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  black,   key="recuperados",  points(tt0,rr0) )$
wxdraw(gk07);

Mf(ii0,ini,fin)*gamma;

gkfh1:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat("Infectados según diferentes vacunas") , 
            point_type = filled_circle, point_size =0.5, 
            color = black, key="sin vacunas", points(tt0,ii0), 
             line_type = dashes, line_width  = 2, 
            color = green,   key="vacuna eficaz",  points(tt0,ve), 
color = red, key="vacuna ineficaz", points(tt0,vi),
color = blue, key="vacuna muy eficaz", points(tt0,vme), yrange = [0, 3500])$
wxdraw(gkfh1);

write_data (makelist([ttl[k],iil0[k],ven[k],vin[k],vmen[k] ],k,1,length(ttl)), "datav3.dat" )$

### Modelo SIRVI

kill(all)$

Ecuaciones del sistema

fs(t,s,i,r, v, i_v) := -alpha*s*(i+i_v)/N+delta*r-nu+mu*v;  /* susceptibles */
fi(t,s,i,r, v, i_v) := alpha*s*(i+i_v)/N-(beta)*i;  /* sanas */
fr(t,s,i,r, v, i_v) := beta*i-delta*r;  /* recuperadas */ 
fv(t,s,i,r, v, i_v) := nu-mu*v-alpha_v*(i+i_v)*v/N+ beta_v*i_v;  /* sanas */
fiv(t,s,i,r, v, i_v) := alpha_v*(i+i_v)*v/N-beta_v*i_v;  /* recuperadas */ 

### Estabilidad

equ1: -alpha*s*(i+i_v)/N+delta*r-nu+mu*v;  /* susceptibles */
equ2: alpha*s*(i+i_v)/N-(beta)*i;  /* sanas */
equ3: beta*i-delta*r;  /* recuperadas */ 
equ4: nu-mu*v-alpha_v*(i+i_v)*v/N+ beta_v*i_v;  /* sanas */
equ5:  alpha_v*(i+i_v)*v/N-beta_v*i_v; /* recuperadas */ 

sc:solve([equ1, equ2, equ3, equ4, equ5], [s,i,r,v,i_v]);

sol1:sc[1];

apply("+",sol1);
ev(%, v+s+r+i_v+i = N);
solve(%, %r1);
s1: ev(sol1, %);

sol2:sc[2];

apply("+",sol2);
ev(%, v+s+r+i_v+i = N);
solve(%, %r2);
s2: ev(sol2, %);

ss2:factor(%);

jac:jacobian([equ1,equ2,equ3,equ4,equ5],[s,i,r, v, i_v]);

jac1:ev(jac,s1),ratsimp;

define(j1(alpha,beta,alpha_v,beta_v,delta,nu,mu,N),jac1);

auto1:eigenvalues(j1(alpha,beta,alpha_v,beta_v,delta,nu,mu,N));

jac2:ev(jac, s2),ratsimp;

define(j2(alpha,beta,alpha_v,beta_v,delta,nu,mu,N),jac2);

define(Gamma1(alpha,beta,alpha_v,beta_v,delta,nu,mu,N), rr/ss );

Gamma1(0.065,0.03,0.04,0.045,0.02,50,0.01,10100);

ratprint: false$
m1:j1(0.065,0.03,0.04,0.045,0.02,50,0.01,10100);
eigenvalues(%),float;

m2:j2(0.065,0.03,0.04,0.045,0.02,50,0.01,10100);
eigenvalues(%),float;

load ("lapack")$
decimal_places(ex,n):=parse_string(printf(false,sconcat("~,",n,"f"),ex));

dgeev (m2)$
eig:%[1]$
for k:1 thru 5 do(print (decimal_places(realpart(eig[k]),16)))$

dgeev (m1)$
eig:%[1]$
for k:1 thru 5 do(print (decimal_places(realpart(eig[k]),16)))$

### SIRVI vacuna eficaz

Mf(lista,p0,pf):= block([ll], ll:makelist(lista[k],k,p0,pf), apply ("+", ll)/length(ll) ) ;
ini:301;fin: 1001;

fs(t,s,i,r, v, i_v) := -alpha*s*(i+i_v)/N+delta*r-nu+mu*v;  /* susceptibles */
fi(t,s,i,r, v, i_v) := alpha*s*(i+i_v)/N-(beta)*i;  /* sanas */
fr(t,s,i,r, v, i_v) := beta*i-delta*r;  /* recuperadas */ 
fv(t,s,i,r, v, i_v) := nu-mu*v-alpha_v*(i+i_v)*v/N+ beta_v*i_v;  /* sanas */
fiv(t,s,i,r, v, i_v) := alpha_v*(i+i_v)*v/N-beta_v*i_v;  /* recuperadas */ 

alpha:0.065$ beta:0.03$ alpha_v : 0.03$ beta_v : 0.06$ delta:0.02 $ nu : 40 $ mu:0.01$
s0:10000$ i0:100$ iv0:0$ v0:0$ N : s0+i0+ v0 + iv0$

realpart(float(eigenvalues(j1(alpha,beta,alpha_v,beta_v,delta,nu,mu,N)))[1]);

aut2:eigenvalues(j2(alpha,beta,alpha_v,beta_v,delta,nu,mu,N)),float$

realpart(aut2[1]);

sol:rk( [fs(t,s,i,r, v, i_v),fi(t,s,i,r, v, i_v),fr(t,s,i,r, v, i_v),fv(t,s,i,r, v, i_v),fiv(t,s,i,r, v, i_v)],   [s,i,r,v,i_v],[s0,i0,0,iv0,v0],[t,0,500,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$
rr0:makelist(sol[k][4],k,1,length(sol),5)$
vv0:makelist(sol[k][5],k,1,length(sol),5)$ 
iv0:makelist(sol[k][6],k,1,length(sol),5)$
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
vvl:append([vac] ,makelist(sol[k][5],k,1,length(sol),50))$
ivl:append([ivac] ,makelist(sol[k][6],k,1,length(sol),50))$
tiv:makelist(sol[k][3],k,1,length(sol),50)+makelist(sol[k][6],k,1,length(sol),50)$
tivn:append([toti],%)$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],vvl[k], ivl[k], tivn[k] ],k,1,length(ttl)), "sirvi1.dat" )$

Mf(ii0+iv0,ini,fin);

gvace:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( "Gráfica con vacunas") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados totales",  points(tt0,ii0+iv0) ,  
            color =  purple,   key="vacunados",  points(tt0,vv0),
 color =  black,   key="recuperados",  points(tt0,rr0))$
wxdraw(gvace);

gvace:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( "Gráfica con vacunas") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  purple,   key="vacunados",  points(tt0,vv0),
    color =  green,   key="infectados vac",  points(tt0,iv0),
color =  black,   key="recuperados",  points(tt0,rr0))$
wxdraw(gvace);

gvaceinf:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( "Gráfica con vacunas") , 
            point_type = filled_circle, point_size =0.5, 
                      line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  green,   key="infectados vac",  points(tt0,iv0),
color =  black,   key="infectados totales",  points(tt0,ii0+iv0))$
wxdraw(gvaceinf);

draw(gvace)$
draw_file(dimensions = [1000,700],  terminal = 'png, file_name = "Grafvac23")$

### SIRVI sin vacuna

alpha:0.065$ beta:0.03$ alpha_v : 0.0$ beta_v : 0.0$ delta:0.02 $ nu : 0 $ mu:0.0$
s0:10000$ i0:100$ iv0:0$ v0:0$ N : s0+i0+ v0 + iv0$

sol:rk( [fs(t,s,i,r, v, i_v),fi(t,s,i,r, v, i_v),fr(t,s,i,r, v, i_v),fv(t,s,i,r, v, i_v),fiv(t,s,i,r, v, i_v)],   [s,i,r,v,i_v],[s0,i0,0,iv0,v0],[t,0,500,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$
rr0:makelist(sol[k][4],k,1,length(sol),5)$
vv0:makelist(sol[k][5],k,1,length(sol),5)$ 
iv0:makelist(sol[k][6],k,1,length(sol),5)$
lll:length(tt0)$

length(sol);

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
vvl:append([vac] ,makelist(sol[k][5],k,1,length(sol),50))$
ivl:append([ivac] ,makelist(sol[k][6],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],vvl[k], ivl[k] ],k,1,length(ttl)), "sirvinv.dat" )$

length(ssl);

ttn:makelist(sol[k][1],k,1,length(sol),50)$
tivnv:makelist(sol[k][3],k,1,length(sol),50)+makelist(sol[k][6],k,1,length(sol),50)$
tivnvn:append([totinv],%)$

Mf(ii0+iv0,ini,fin);

novac:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( "Gráfica con vacunas") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  purple,   key="vacunados",  points(tt0,vv0),
    color =  green,   key="infectados vac",  points(tt0,iv0),
color =  black,   key="recuperados",  points(tt0,rr0))$
wxdraw(novac);

compvenv:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat(" Infectados sin vacunas (negro), con vacunas (rojo)") , 
            point_type = filled_circle, point_size =0.5, 
            color = black, key="sin vacunas", points(ttn,tivnv), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="con vacunas",  points(ttn,tiv))$
wxdraw(compvenv);

write_data (makelist([ttl[k],tivnvn[k],tivn[k] ],k,1,length(ttl)), "sirvicomp1.dat" )$

### SIRVI Vacuna poco eficaz

alpha:0.065$ beta:0.03$ alpha_v : 0.03$ beta_v : 0.06$ delta:0.02 $ nu : 40 $ mu:0.1$
s0:10000$ i0:100$ iv0:0$ v0:0$ N : s0+i0+ v0 + iv0$

realpart(float(eigenvalues(j1(alpha,beta,alpha_v,beta_v,delta,nu,mu,N)))[1]);

realpart(float(eigenvalues(j2(alpha,beta,alpha_v,beta_v,delta,nu,mu,N)))[1]);

sol:rk( [fs(t,s,i,r, v, i_v),fi(t,s,i,r, v, i_v),fr(t,s,i,r, v, i_v),fv(t,s,i,r, v, i_v),fiv(t,s,i,r, v, i_v)],   [s,i,r,v,i_v],[s0,i0,0,iv0,v0],[t,0,500,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$
rr0:makelist(sol[k][4],k,1,length(sol),5)$
vv0:makelist(sol[k][5],k,1,length(sol),5)$ 
iv0:makelist(sol[k][6],k,1,length(sol),5)$
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
vvl:append([vac] ,makelist(sol[k][5],k,1,length(sol),50))$
ivl:append([ivac] ,makelist(sol[k][6],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],vvl[k], ivl[k] ],k,1,length(ttl)), "sirvipe.dat" )$

tivpe:makelist(sol[k][3],k,1,length(sol),50)+makelist(sol[k][6],k,1,length(sol),50)$
tivnpe:append([totipe],%)$

Mf(ii0+iv0,ini,fin);

vacpe:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( "Gráfica con vacunas") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  purple,   key="vacunados",  points(tt0,vv0),
    color =  green,   key="infectados vac",  points(tt0,iv0),
color =  black,   key="recuperados",  points(tt0,rr0))$
wxdraw(vacpe);

### SIRVI Ritmo lento de Vacunación

alpha:0.065$ beta:0.03$ alpha_v : 0.03$ beta_v : 0.06$ delta:0.02 $ nu : 5 $ mu:0.01$
s0:10000$ i0:100$ iv0:0$ v0:0$ N : s0+i0+ v0 + iv0$

realpart(float(eigenvalues(j1(alpha,beta,alpha_v,beta_v,delta,nu,mu,N)))[1]);

realpart(float(eigenvalues(j2(alpha,beta,alpha_v,beta_v,delta,nu,mu,N)))[1]);

sol:rk( [fs(t,s,i,r, v, i_v),fi(t,s,i,r, v, i_v),fr(t,s,i,r, v, i_v),fv(t,s,i,r, v, i_v),fiv(t,s,i,r, v, i_v)],   [s,i,r,v,i_v],[s0,i0,0,iv0,v0],[t,0,500,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$
rr0:makelist(sol[k][4],k,1,length(sol),5)$
vv0:makelist(sol[k][5],k,1,length(sol),5)$ 
iv0:makelist(sol[k][6],k,1,length(sol),5)$
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
vvl:append([vac] ,makelist(sol[k][5],k,1,length(sol),50))$
ivl:append([ivac] ,makelist(sol[k][6],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],vvl[k], ivl[k] ],k,1,length(ttl)), "sirvirl.dat" )$

tivrl:makelist(sol[k][3],k,1,length(sol),50)+makelist(sol[k][6],k,1,length(sol),50)$
tivnrl:append([totirl],%)$

Mf(ii0+iv0,ini,fin);

vacrl:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( "Gráfica con vacunas") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  purple,   key="vacunados",  points(tt0,vv0),
    color =  green,   key="infectados vac",  points(tt0,iv0),
color =  black,   key="recuperados",  points(tt0,rr0))$
wxdraw(vacrl);

### SIRVI Vacuna muy eficaz

alpha:0.065$ beta:0.03$ alpha_v : 0.03$ beta_v : 0.06$ delta:0.02 $ nu : 40 $ mu:0.005$
s0:10000$ i0:100$ iv0:0$ v0:0$ N : s0+i0+ v0 + iv0$

realpart(float(eigenvalues(j1(alpha,beta,alpha_v,beta_v,delta,nu,mu,N)))[1]);

realpart(float(eigenvalues(j2(alpha,beta,alpha_v,beta_v,delta,nu,mu,N)))[1]),float;

sol:rk( [fs(t,s,i,r, v, i_v),fi(t,s,i,r, v, i_v),fr(t,s,i,r, v, i_v),fv(t,s,i,r, v, i_v),fiv(t,s,i,r, v, i_v)],   [s,i,r,v,i_v],[s0,i0,0,iv0,v0],[t,0,500,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$
rr0:makelist(sol[k][4],k,1,length(sol),5)$
vv0:makelist(sol[k][5],k,1,length(sol),5)$ 
iv0:makelist(sol[k][6],k,1,length(sol),5)$
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
vvl:append([vac] ,makelist(sol[k][5],k,1,length(sol),50))$
ivl:append([ivac] ,makelist(sol[k][6],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],vvl[k], ivl[k] ],k,1,length(ttl)), "sirvime.dat" )$

tivme:makelist(sol[k][3],k,1,length(sol),50)+makelist(sol[k][6],k,1,length(sol),50)$
tivnme:append([totime],%)$

Mf(ii0+iv0,ini,fin);

vacme:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( "Gráfica con vacunas") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  purple,   key="vacunados",  points(tt0,vv0),
    color =  green,   key="infectados vac",  points(tt0,iv0),
color =  black,   key="recuperados",  points(tt0,rr0))$
wxdraw(vacme);

### SIRVI Ritmo alto de Vacunación

alpha:0.065$ beta:0.03$ alpha_v : 0.03$ beta_v : 0.06$ delta:0.02 $ nu : 75 $ mu:0.01$
s0:10000$ i0:100$ iv0:0$ v0:0$ N : s0+i0+ v0 + iv0$

realpart(float(eigenvalues(j1(alpha,beta,alpha_v,beta_v,delta,nu,mu,N)))[1]);

realpart(float(eigenvalues(j2(alpha,beta,alpha_v,beta_v,delta,nu,mu,N)))[1]),float;

sol:rk( [fs(t,s,i,r, v, i_v),fi(t,s,i,r, v, i_v),fr(t,s,i,r, v, i_v),fv(t,s,i,r, v, i_v),fiv(t,s,i,r, v, i_v)],   [s,i,r,v,i_v],[s0,i0,0,iv0,v0],[t,0,500,0.1] )$
ll:length(sol);
tt0:makelist(sol[k][1],k,1,length(sol),5)$ 
ss0:makelist(sol[k][2],k,1,length(sol),5)$  
ii0:makelist(sol[k][3],k,1,length(sol),5)$
rr0:makelist(sol[k][4],k,1,length(sol),5)$
vv0:makelist(sol[k][5],k,1,length(sol),5)$ 
iv0:makelist(sol[k][6],k,1,length(sol),5)$
lll:length(tt0)$

/* Para exportar los datos y pintarlos con TikZ en LaTeX */
ttl:append([time] ,makelist(sol[k][1],k,1,length(sol),50))$
ssl:append([suc] ,makelist(sol[k][2],k,1,length(sol),50))$
iil:append([inf] ,makelist(sol[k][3],k,1,length(sol),50))$
rrl:append([rec] ,makelist(sol[k][4],k,1,length(sol),50))$
vvl:append([vac] ,makelist(sol[k][5],k,1,length(sol),50))$
ivl:append([ivac] ,makelist(sol[k][6],k,1,length(sol),50))$
write_data (makelist([ttl[k],ssl[k],iil[k],rrl[k],vvl[k], ivl[k] ],k,1,length(ttl)), "sirvira.dat" )$

tivra:makelist(sol[k][3],k,1,length(sol),50)+makelist(sol[k][6],k,1,length(sol),50)$
tivnra:append([totira],%)$

Mf(ii0+iv0,ini,fin);

vacra:gr2d( grid = true, line_type =short_long_dashes, line_width  = 2, nticks=100,
        xaxis =true, xaxis_type=solid, yaxis =true, yaxis_type=solid, 
            xlabel="unidades de tiempo", 
            title=concat( "Gráfica con vacunas") , 
            point_type = filled_circle, point_size =0.5, 
            color = blue, key="susceptibles", points(tt0,ss0), 
             line_type = dashes, line_width  = 2, 
            color = red,   key="infectados",  points(tt0,ii0) ,  
            color =  purple,   key="vacunados",  points(tt0,vv0),
    color =  green,   key="infectados vac",  points(tt0,iv0),
color =  black,   key="recuperados",  points(tt0,rr0))$
wxdraw(vacra);

write_data (makelist([ttl[k],tivnvn[k],tivn[k], tivnpe[k],tivnrl[k],tivnme[k],tivnra[k] ],k,1,length(ttl)), "sirvicomp.dat" )$