%% Assignment 2
%
% (c) Konstantin Borschiov, 2008
% Estimate population growth using Numerical and Symbolic methods
%

% Prepare the workspace
clear all;
clc;
clf;

%% Population when b = 0

dt = 0.001;
start = 0;
stop = 10;
time = [start:dt:stop];

% Initialize vector for storing numerical solution
N1 = zeros(size(time));
N1(1) = 10;

% Initialize population growth modifiers
a = 0.5;
b = 0;

% Numerical solution
for i = 2:length(N1)
    N1(i) = N1(i-1) + a*N1(i-1)*dt;
end

% Initialize vector for storing symbolic solution
N2 = zeros(size(N1));
N2(1) = N1(1);

% Symbolic solution
N2 = N2(1).*exp(a*time);


% Plot Numeric and Symbolic solution on a graph
figure(1);
subplot(211)
plot (time, N1,'*b', time, N2, '+g');
grid on;
xlabel('Time');
ylabel('Population');
title(['Population growth with startin gpopulation of ' num2str(N1(1)) ', a = ' num2str(a) ' , b= ' num2str(b) ' ']);
legend('Numerical Solution', 'Symbolic Solution', 'Location', 'Best');

% Plot difference between numeric and symbolic solution
subplot(212);
plot (time, N2-N1);
xlabel('Time');
ylabel('Difference');
title('Difference between Numerical and Symbolic solution');



%% Population when b ~= 0


% Initialize vector for storing numerical solution
N3 = zeros(size(time));
N3(1) = 1000;

% Initialize population growth modifiers
a2 = 15;
b2 = 0.01;

% Numerical Solution
for i = 2:length(N3)
    temp1 = a2*N3(i-1);
    temp2 = b2*(N3(i-1))^2;
    N3(i) = N3(i-1) + (temp1 - temp2)*dt;
end

% Initialize vector for storing symbolic solution
N4 = zeros(size(N3));
N4(1) = N3(1);
temp3 = zeros(size(N3)); % Temporary array to simplify expression

% Symbolic solution
temp3 = (N4(1)*exp(a2*time))./(a2 - b2*N4(1)); % Temporary array to simplify expression
N4 = temp3*a2./(1+temp3*b2);


% Plot Numeric and Symbolic solution on a graph
figure(2);
subplot(211);
plot(time, N3, '-b', time, N4, '-.r');
grid on;
xlabel('Time');
ylabel('Population');
title(['Population growth with starting population of ' num2str(N3(1)) ', a = ' num2str(a2) ' , b= ' num2str(b2) ' ']);
legend('Numerical Solution', 'Symbolic Solution', 'Location', 'Best');

% Plot difference between numeric and symbolic solution
subplot(212);
plot (time, N4-N3);
xlabel('Time');
ylabel('Difference');
title('Difference between Numerical and Symbolic solution');