sinan

The generated figures are added.

clear all; clc; close all; clf;
sinan = sinan();

t0 = 0; tfinal = 60; N = 100;
h = (tfinal-t0)/N;
t = t0:h:tfinal;
%% Initial conditions
y0=[100, 10, 20];
y01=[100];

%% Parameters
Lambda = 0.1; beta = 0.4; mu = 0.01; gamma = 0.7; % parameters for model 1
a = 0.2; % parameter for model 2

Models

SIR Epidemic model

subject to

Decay model

subject to

%% Models
model_1 = @(t, X) [Lambda- beta*X(1)*X(2);
    beta*X(1)*X(2)-(gamma+mu)*X(2);
    mu*X(2)];
model_2 = @(t, y) -a*y;

% fractional order
alpha = 1;

% Choose a solver
solver_1 = 'caputo';
solver_2 = 'abc';
solver_3 = 'CLASSICAL';

% How to solve a system of equations?
[t_euler_caputo_1, y_euler_caputo_1] = sinan.euler(model_1, alpha, t, h, y0, solver_1);
[t_euler_ABC_1, y_euler_ABC_1] = sinan.euler(model_1, alpha, t, h, y0, solver_2);
[t_euler_clas_1, y_euler_clas_1] = sinan.euler(model_1, alpha, t, h, y0, solver_3);
[t_rk4_caputo_1, y_rk4_caputo_1] = sinan.rk4(model_1, alpha, t, h, y0, solver_1);
[t_rk4_ABC_1, y_rk4_ABC_1] = sinan.rk4(model_1, alpha, t, h, y0, solver_2);
[t_rk4_clas_1, y_rk4_clas_1] = sinan.rk4(model_1, alpha, t, h, y0, solver_3);

% How to solve an equation?
[t_euler_caputo_2, y_euler_caputo_2] = sinan.euler(model_2, alpha, t, h, y01, solver_1);
[t_euler_ABC_2, y_euler_ABC_2] = sinan.euler(model_2, alpha, t, h, y01, solver_2);
[t_euler_clas_2, y_euler_clas_2] = sinan.euler(model_2, alpha, t, h, y01, solver_3);
[t_rk4_caputo_2, y_rk4_caputo_2] = sinan.rk4(model_2, alpha, t, h, y01, solver_1);
[t_rk4_ABC_2, y_rk4_ABC_2] = sinan.rk4(model_2, alpha, t, h, y01, solver_2);
[t_rk4_clas_2, y_rk4_clas_2] = sinan.rk4(model_2, alpha, t, h, y01, solver_3);

figure(1)
plot(t_euler_caputo_1, y_euler_caputo_1(:, 1), 'r', t_euler_ABC_1, y_euler_ABC_1(:, 1), 'b', t_euler_clas_1, y_euler_clas_1(:, 1), 'k')
legend('Euler caputo', 'Euler abc', 'Euler classical')

figure(2)
plot(t_rk4_caputo_1, y_rk4_caputo_1, 'r', t_rk4_ABC_1, y_rk4_ABC_1, 'b', t_rk4_clas_1, y_rk4_clas_1, 'k')
legend('rk4 caputo', 'rk4 abc', 'rk4 classical')
figure(3)
plot(t_euler_clas_1, y_euler_clas_1, 'r', t_rk4_clas_1, y_rk4_clas_1, 'b')
legend('euler classical', 'rk4 classical')

figure(4)
plot(t_euler_caputo_2, y_euler_caputo_2, 'r', t_euler_ABC_2, y_euler_ABC_2, 'b', t_euler_clas_2, y_euler_clas_2, 'k')
legend('Euler caputo', 'Euler abc', 'Euler classical')
xlabel('time (t)')
ylabel('y(t)')

figure(5)
plot(t_rk4_caputo_2, y_rk4_caputo_2, 'r', t_rk4_ABC_2, y_rk4_ABC_2, 'b', t_rk4_clas_2, y_rk4_clas_2, 'k')
legend('rk4 caputo', 'rk4 abc', 'rk4 classical')
xlabel('time (t)')
ylabel('y(t)')

figure(6)
plot(t_euler_clas_2, y_euler_clas_2, 'r', t_rk4_clas_2, y_rk4_clas_2, 'b')
legend('euler classical', 'rk4 classical')
xlabel('time (t)')
ylabel('y(t)')