According to the PDF, Eq. (6) is the yield distribution function 
. Though I'm unfamiliar with this model, I'd advise you to type out all equations and paramters (constants) in MATLAB code in the odefcn() function model, like the following Van de Pol example. This allows other people to conveniently test the simulation and correct your code. It's okay even if it does not work.
. Though I'm unfamiliar with this model, I'd advise you to type out all equations and paramters (constants) in MATLAB code in the odefcn() function model, like the following Van de Pol example. This allows other people to conveniently test the simulation and correct your code. It's okay even if it does not work.[t, y] = ode45(@odefcn, [0 20], [2; 0]);
plot(y(:,1), y(:,2)), grid on
title('Limit cycle of Van der Pol oscillator (\mu = 1)');
xlabel('y_{1}'), ylabel('y_{2}')
% Van der Pol oscillator
function dydt = odefcn(t, y)
% parameter
mu = 1;
% ODE in state-space form
dydt = [y(2); % ode dy1/dt for state y1
mu*(1 - y(1)^2)*y(2) - y(1)]; % ode dy2/dt for state y2
end
