function yourIntegration
tResult = [];
xResult = [];
tStep = [7 14 21 28 35 42 49 56 63 70 77 84];
x0 = [64700 0 0.0033];
for index = 2:numel(tStep)
beta = 0.43e-08 + (4.28e-08 - 0.43e-08) * exp(-0.20*t(index - 1))
af = @(t,x) f(t, x, beta);
t = tStep(index-1:index);
[t, x] = ode45(af, t, x0);
tResult = cat(1, tResult, t);
xResult = cat(1, xResult, x);
x0 = x(end, :);
end
function dx = f(t,x, beta)
dx = [3494-0.054*x(1)-beta*x(1)*x(3); ...
beta*x(1)*x(3) - 0.41*x(2); ...
50000*x(2) - 23*x(3)];