Matlab:Basic population model

%
% BASIC MAGICAL GIRL-WITCHES MODEL
% This file has been produced for Matlab 7.4.0
% It is 100% compatible with FreeMat 4.0, however, you may want to change the plot commands to
% plot(M, 'g.-') and plot(W, 'r.-') for aesthetics purpose
%
% If you want to modify the constants, please bear these restrictions in mind:
% -C should be an integer
% -D, B, F and K should be values in [0,1]
% -B should be greater than K (or else, the number of witches will be negative)
%


C = 100;     % Amount of Magical Girls contracted by Kyubey, per iteration
D = 0.25;    % Proportion of girls getting killed, per iteration
B = 0.125;   % Proportion of girls becoming witches, per iteration
F = 0.015;   % Proportion of familiars maturing into witches, per iteration
K = 0.12;    % Proportion of witches getting killed

steps = 100; % Number of iterations for the simulation

M(1) = 0;    % Number of magical girls at first
W(1) = 0;    % Number of witches at first

hold on

for t=1:steps
    deltaM = C-(D+B)*M(t);
    deltaW = F*W(t) + (B-K)*M(t);
    
    M(t+1)=M(t)+deltaM;
    W(t+1)=W(t)+deltaW;
end

title('Evolution of the Magical Girls population over time')
xlabel('time')
ylabel('number of magical girls/witches')
plot(M, 'g-')
plot(W, 'r-')
legend('Magical Girls', 'Witches', 'location', 'southeast')
legend('boxoff')

hold off