clear; format long; global N C0 ts;
p = 8; ts = [1:35]; given = ts(end);
Tdata = ts(given+1-p:given); Tdata = Tdata-Tdata(1);
Cs =[278 309 571 830 1297 1985 2741 4537 5997 7736 9720 11821 14411 17238 20471 ...
24363 28060 31211 34598 37251 40235 42708 44730 46550 48548 50054 51174 ...
70635 72528 74280 74675 75569 76392 77042 77262];
Cdata = Cs(given+1-p:given); p = Tdata(end);
N = 1.0e+8; C(1) = Cdata(1); C0 = C(1); b0 = 1.0; g0 = 1.0; Us = 0.1*C0;
param = lsqcurvefit(@SUCmodel,[b0 g0 Us],Tdata,Cdata,...
[1.0e-3 1.0e-3 0.01*C0],[10 10 5*C0]);
b = param(1); g = param(2); U(1) = param(3); S(1) = N-U(1)-C(1);
Nt = 1000*p; dt = 0.001;
for i = 1:Nt
S(i+1) = S(i) - dt*b*S(i)*U(i)/N;
U(i+1) = U(i) + dt*(b*S(i)*U(i)/N - g*U(i));
end
C = N-S-U;
figure(1); clf; grid on; hold on
t=linspace(0,p,Nt+1);
plot(Tdata,Cdata,'ko','linewidth',2);
plot(t,C,'k','linewidth',1.5); plot(t,U,'k:','linewidth',2);
legend('Confirmed case data','C - Best fitting confirmed', ...
'U - Unidentified infected','Location','best')
axis([t(1) t(end) 0 1.1*max(U(end),C(end))]); xlabel('t'); ylabel('population')


function f = SUCmodel(Parameter, Tdata)
global N C0;
b = Parameter(1);
g = Parameter(2);
U(1) = Parameter(3);
C(1) = C0;
S(1) = N-C(1)-U(1);
Nt = 1000*Tdata(end);
dt = 0.001;
for i = 1:Nt
S(i+1) = S(i) - dt*b*S(i)*U(i)/N;
U(i+1) = U(i) + dt*(b*S(i)*U(i)/N - g*U(i));
end
C = N-S-U;
t = linspace(0,Tdata(end),Nt+1);
f = interp1(t,C,Tdata);
end