Are gamma and beta parameters that result from the gamma / beta functions? Or just scalar parameters? I suggest using other names if they are scalars, because Matlab inbuilt functions are called like that, which will produce errors.
Run this part only once - ideally in a seperate script.Save or copy the result into another script and optimize then without the symbolic calculation. Therefore use beta (and maybe the others) as an additional input.
syms S(t) I(t) R(t) D(t) Beta pop Gamma mu
odeS = diff(S) == -Beta*I*S/pop;
odeI = diff(I) == Beta*I*S/pop - Gamma*I-mu*I;
odeR = diff(R) == Gamma*I;
odeD = diff(D) == mu*I;
% Transform the equations for the numerical solver
odes = [odeS; odeI; odeR; odeD];
odes2 = odeToVectorField(odes);
eq_mat = matlabFunction(odes2, 'Vars', {'t', 'Y', 'Beta', 'pop', 'Gamma', 'mu'})
The numerical solution process is fast, but symbolic calculations are not. So this part should be the only one that runs repeated durng the optimization:
eq_mat = @(t,Y,Beta,pop,Gamma,mu)[-Gamma.*Y(1)-mu.*Y(1)+(Beta.*Y(1).*Y(2))./pop;-(Beta.*Y(1).*Y(2))./pop;Gamma.*Y(1);mu.*Y(1)]
s0 = 1;
i0 = -1;
r0 = 1;
d0 = -1;
ic = [s0, i0, r0, d0];
Gamma = 1;
pop = 2;
mu = -1;
Beta = 0.5;
tspan = [0, 100];
[t, y]= ode45(@(t,Y)eq_mat(t,Y,Beta,pop,Gamma,mu), tspan, ic);
plot(t,y)
