Error using nlinfit (Check FunVals)

Hi there.. i am trying to fit my odel equations to some data but i keep on getiing error.. please help

below is my code

% Data fitting for the covid-19 outbreak in Lagos

covid19_nlinfit

 
                                     Norm of         Norm of
   Iteration             SSE        Gradient           Step 
  -----------------------------------------------------------
           0     8.39013e+06

Warning: Rank deficient, rank = 2, tol = 1.269216e-11.

           1     7.37178e+06     5.39372e+06         2.01582

Warning: Failure at t=1.591672e+01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.842171e-14) at time t.

Error using nlinfit>checkFunVals (line 649)
The function you provided as the MODELFUN input has returned Inf or NaN values.

Error in nlinfit>LMfit (line 596)
if funValCheck && ~isfinite(sse), checkFunVals(r); end

Error in nlinfit (line 284)
[beta,J,~,cause,fullr] = LMfit(X,yw, modelw,beta,options,verbose,maxiter);

Error in solution>covid19_nlinfit (line 36)
[values] = nlinfit(day, Totalcases, @covidinfection, [a(1) a(2) a(3)], options);

function covid19_nlinfit
clear all
clc

global a x1 beta_c ep_g c_g sigma nu theta psi d_o d_d gamma_i alpha ...
    delta ep_s ep_i ep_a E0 A0 I0 

ep_g = 0.9; c_g = 0.9; delta = 0; ep = 0; ep_s = ep; ep_a=ep; ep_i=ep;
sigma = 1/5.2; nu = 0.5; 
gamma_i = 1/15; 
gamma_a = 0.13978; gamma_o = 0.13978;
d_o = 0.015; d_d = 0.015;
alpha = 0.5; 

%==========================================================================
% Lagos data, starting from March 16, 2020 till May 3, 2020.
day = [1 3 4 6 7 8 9 10 11 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 ...
    43 44 45 46 47 48];

Totalcases = [1 5 9 19 22 25 29 32 44 59 68 81 82 91 98 109 109 120 120 130 145 158 166 174 176 189 214 232 251 283 ...
    306 376 376 430 504 582 657 689 731 764 844 931 976 1006 1068];
Activecases = [1 5 8 19 22 25 29 32 43 58 67 75 76 85 81 86 86 91 87 96 111 117 118 119 116 123 139 140 154 182 200 ...
    266 266 309 383 460 525 548 584 602 682 718 756 753 791];
%==========================================================================
%Initial guesses for parameters to be estimated
beta_c=0; theta=0;psi=0;
%    E0 = 16; A0 = 20; I0 = 10;

%======================================================================
% Using NLINFIT function
a = [beta_c theta psi]; % E0 A0 I0;
x1 = day;
options = statset('Display', 'iter');
[values] = nlinfit(day, Totalcases, @covidinfection, [a(1) a(2) a(3)], options);
%[values] = nlinfit(day, Activecases, @covidinfection, [a(1) a(2) a(3) a(4) a(5) a(6)], options);   
beta_c = values(1)
theta  = values(2)
psi    = values(3)
%     E0 = values(4)
%     A0 = values(5)
%     I0 = values(6)
%=====================================================================
%Control Reproduction Number
R_c = beta_c*(1-ep_g*c_g)*(1-delta)*(1-ep_s)*((alpha*nu*(1-ep_a)/...
    (theta+gamma_a)) + ((1-nu)*(1-ep_i)/(psi+d_o+gamma_o))); 

[t,y] = ode15s(@covidmodel1, [0,50],[14368332,16,20,10,1,0,0]); %E0 A0 I0;

plot(t,y(:,7),'-b',day,Totalcases,'rs')
xlabel('Time (days)'),ylabel('Cumulative number of reported cases')
legend('Model','Lagos COVID-19 Data')
%plot(t,y(:,5),'-b',day,Activecases,'rs')

%========================================================================
%  Using NLINFIT function for data fitting
    function y = covidinfection(a,x1)  
        
        beta_c = a(1);
        theta = a(2);
        psi   = a(3);
        %      E0 = a(4);
        %      A0 = a(5);
        %      I0 = a(6);
        
        [t,y] = ode15s(@covidmodel1, [0,x1(end)],[14368332,16,20,10,1,0,0]); %E0 A0 I0;
        
        y = interp1(t,y(:,7),x1); % cumulative reported cases
        %y = interp1(t,y(:,5),x1); % active cases
    end
%=========================================================================

    function dx = covidmodel1(t,x)
        
        dx = [0;0;0;0;0;0;0];
        N = x(1)+x(2)+x(3)+x(4)+x(5)+x(6);
        beta = beta_c*(1-ep_g*c_g)*(1-delta)*(1-ep_s);
        
        dx(1) = -beta*x(1)*(alpha*(1-ep_a)*x(3)+(1-ep_i)*x(4))/N; 
        dx(2) = beta*x(1)*(alpha*(1-ep_a)*x(3)+(1-ep_i)*x(4))/N - sigma*x(2);
        dx(3) = nu*sigma*x(2) - theta*x(3) - gamma_a*x(3);
        dx(4) = (1-nu)*sigma*x(2) - psi*x(4) - d_o*x(4) - gamma_o*x(4);
        dx(5) = theta*x(3) + psi*x(4) - gamma_i*x(5) -d_d*x(5); % Active cases, detected (asymptomatic and symptomatic)
        dx(6) = gamma_i*x(5) + gamma_a*x(3) + gamma_o*x(4);
        dx(7) = theta*x(3) + psi*x(4); % Cumulative number of new cases reported (asymptomatic and symptomatic)
    end

end

Answer 1

[Matt J]

For one thing,

y = interp1(t,y(:,7),x1,'spline','extrap'); % cumulative reported cases

For another, some sort of bounds or constraints on the variables are needed to prevent the ODE solver from considering ill-behaved systmes of equations. nlinfit does not let you impose such constraints, so you may have to look to another solver, e.g. lsqcurvefit.