MATLAB: SIRD model- Input arguments errors!

argumentsinput

I am trying to calculate a solution to a SIRD model using an ERK method but I keep getting a number of different errors no matter what I try to change. Can anyone see any way in which this programme will run and give me an answer? I have tried adding all the parameters to the stage values k1-4 but the programme just bombs and won't even come back with an error it will just never stop running.
This is my code:
clear
clc
% Define ODE function
f = @(t, y, N, beta, gamma, mu) [ -(beta*y(2)*y(1))/N, (beta*y(2)*y(1))/N-gamma*y(2)-mu*y(2), gamma*y(2), mu*y(2) ];
% Define IVP
tspan = [ 0, 3 ];
y0 = [ 4990000, 10000, 0, 0 ];
h = 0.5;
N = 5000000;
beta = 0.4;
gamma = 1/21;
mu = 2/1000;
% ERK
y(1, :) = y0;
t(1)=tspan(1);
nsteps = (tspan(2)-tspan(1)) / h;
for n = 1 : nsteps
    k1 = f(t(n), y(n, :));
    k2 = f(t(n) + h/2, y(n, :) + h*k1/2);
    k3 = f(t(n) + 3*h/4, y(n, :) + 3*h*k2/4);
    k4 = f(t(n) + h, y(n, :) + h*k3);
    y(n+1, :) = y(n, :) + (h/36)*(7*k1+18*k2+8*k3+3*k4);
    t(n+1) = t(n) + h;
end
And this is my error:
Not enough input arguments.
Error in Q3b>@(t,y,N,beta,gamma,mu)[-(beta*y(2)*y(1))/N,(beta*y(2)*y(1))/N-gamma*y(2)-mu*y(2),gamma*y(2),mu*y(2)] (line
5)
f = @(t, y, N, beta, gamma, mu) [ -(beta*y(2)*y(1))/N, (beta*y(2)*y(1))/N-gamma*y(2)-mu*y(2), gamma*y(2), mu*y(2) ];
Error in Q3b (line 21)
k1 = f(t(n), y(n, :));

Best Answer

f = @(t, y, N, beta, gamma, mu) [ -(beta*y(2)*y(1))/N, (beta*y(2)*y(1))/N-gamma*y(2)-mu*y(2), gamma*y(2), mu*y(2) ];
That defines an anonymous function that expects 6 input parameters.
tspan = [ 0, 3 ];
y0 = [ 4990000, 10000, 0, 0 ];
h = 0.5;
N = 5000000;
beta = 0.4;
gamma = 1/21;
mu = 2/1000;
That defines some constants, but does not substitute the constants into the function handle. For everything except display purposes (and checking to see if you passed too many inputs), your function handle
f = @(t, y, N, beta, gamma, mu) [ -(beta*y(2)*y(1))/N, (beta*y(2)*y(1))/N-gamma*y(2)-mu*y(2), gamma*y(2), mu*y(2) ];
is equivalent to
f = @(varargin) [ -(varargin{4}*varargin{2}(2)*varargin{2}(1))/varargin{3}, (varargin{4}*varargin{2}(2)*varargin{2}(1))/varargin{3}-varargin{5}*varargin{2}(2)-varargin{6}*varargin{2}(2), varargin{5}*varargin{2}(2), varargin{6}*varargin{2}(2) ];
nothing you assign to any variable mentioned in the function handle definition makes any definition to the execution of the function handle.
Furthermore, your beta, gamma, and so on, are named parameters in the function call definition, and named parameters always override captured variables:
gamma = 123;
delta = 456;
F = @(alpha, beta, gamma) alpha+beta+gamma+delta
that function handle definition "captures" the value of delta so that any changes to delta would not affect the value bundled into the function, but gamma here in the function body would be hard-bound to the position parameters, like
F is defined as a function with three parameters whose body is (first positional parameter passed) + (second positional parameter passed) + (third positional parameter passed) + (saved copy of delta)
and that is a problem for you, as you are only passing in two parameters to f in your lines such as
    k1 = f(t(n), y(n, :));
You have two choices:
  1. You can change how you call f to always pass in those parameters, such as |k1 = f(t(n), y(n,:), N, beta, gamma, nu); OR
  2. You can assign to those variables first before you define the function handle, and do not name the constants as parameters: @(t,y) [ -(beta*y(2)*y(1))/N, etc counting on the anonymous handle to "capture" the values of the variables.
Related Question