function HW1_solutions
%
% Solutions to Homework 1, 2014
% 10 POINTS + 10 extra credit
%
    close all
%
%%   Problem 9 [1 POINT]
%
    v = randi(100,1000,1);
    w = randi(100,1000,1);
 %
 %% Problem 10 [3 POINTS]
 %
    n_equal_vals_in_v_and_w = count_equal_vals(v,w);
    fprintf('%s, %d\n','number of equal values in v and w:',n_equal_vals_in_v_and_w)
 %
 %% Problem 11 [3 POINTS + 2 extra credit for labeling graphs and/or eliminating the error]
 %
    a = 10;
    k = -1;
    r0 = 1;
    tstop = 40;
    [time_vals,sol_vals] = ode45(@universe_equation,[0,tstop],r0);
    figure11 = figure;
    axes11 = axes('Parent',figure11,'FontSize',14);
    box(axes11,'on');
    hold(axes11,'all');
    plot(time_vals,sol_vals,'Parent',axes11,'LineWidth',2,'Color',[1 0 0],...
        'DisplayName','k=-1 (inflationary universe)');
    hold on
% Solution with k=+1, default tolerance.   Note that the solution
% is in correct - for this case the solution should asymptote to
% r=a/k but it does not (unless the ode equation is modified to enforce
% this
    k=1;
    [time_vals,sol_vals] = ode45(@universe_equation,[0,tstop],r0);
    plot(time_vals,sol_vals,'Parent',axes11,'LineWidth',2,...
    'DisplayName','k=1, default tolerance - incorrect)');
%   Solution with k=1 and reduced tolerance.   This gives a better
%   solution although there are still some glitches (see the error message)
    options = odeset('RelTol',0.000001);
    [time_vals,sol_vals] = ode45(@universe_equation,[0,tstop],r0,options);
    plot(time_vals,sol_vals,'Parent',axes11,'LineWidth',2,'Color',[0 1 0],...
    'DisplayName','k=1, reduced tolerance- better)');
    xlabel('Time (au)','FontSize',14);
    ylabel('Size of universe (au)','FontSize',14);
    legend1 = legend(axes11,'show');
    set(legend1,...
    'Position',[0.20581854043393 0.639987466926612 0.395266272189349 0.235528942115768]);
    
    function drdt = universe_equation(t,r)
 %      Uncomment the conditional statement to stop the peculiar
 %      error message
 %      if (a/r-k>0)
        drdt = sqrt(a/r-k);
 %      else
 %      drdt = 0;
 %     end
    end
    
%%
% Problem 12    [3 POINTS + 2 extra credit for labeling graphs]
    beta = 0.4;
    gamma = 0.3;
    S0 = 0.99;
    I0 = 0.01;
    R0 = 0;
    tstop = 90;
    initial_w = [S0;I0;R0];
    options = odeset('RelTol',0.000001);
    [time_vals,sol_vals] = ode45(@SIR_equation,[0,tstop],initial_w,options);
    figure1 = figure;
    axes1 = axes('Parent',figure1,'FontSize',14);
    box(axes1,'on');
    hold(axes1,'all');
    plot1 = plot(time_vals,sol_vals,'Parent',axes1,'LineWidth',2);
    set(plot1(1),'DisplayName','S');
    set(plot1(2),'DisplayName','I');
    set(plot1(3),'DisplayName','R');
    xlabel('Time (days)','FontSize',14);
    ylabel('Population Fraction','FontSize',14);
    legend1 = legend(axes1,'show');
    set(legend1,...
    'Position',[0.743495336278841 0.755798542080848 0.154639175257732 0.156394963552021]);


    function dwdt = SIR_equation(t,w)
        S = w(1);  I = w(2);
        dSdt = -beta*S*I;
        dIdt = beta*S*I-gamma*I;
        dRdt = gamma*I;
        dwdt = [dSdt;dIdt;dRdt];
    end

%% Problem 13 [6 POINTS extra credit]

n = 1000;
mygreatlistoprimes = findprimes(n);

fprintf('%s %d \n ','Here is a list of primes up to ',n);
fprintf('%d %d %d %d %d %d %d %d %d %d \n ',mygreatlistoprimes);
fprintf('\n')


end
function n_equal_vals = count_equal_vals(vector1,vector2)
% Function to count number of equal values in two vectors
n_vals = min(length(vector1),length(vector2));
n_equal_vals = 0;

for i = 1:n_vals
    if (vector1(i)==vector2(i))
        n_equal_vals = n_equal_vals + 1;
    end
end




end
function primes = findprimes(maxnumber)
% Function to calculate all the primes less than maxnumber

% The vector allnums will have the following properties
% allnums(i) = 0 if i is prime
% allnums(i) = n where n is the lowest prime factor of i otherwise

allnums = zeros(maxnumber+1,1);  % Assume everything is prime
currentprime = 2;  % Start with the lowest prime number
while(currentprime<maxnumber+1) % Keep looping until the current prime exceeds maxnumber
    notprime = 2*currentprime;
    count = 2;
    while (notprime<maxnumber+1) % Mark all the numbers that are factors of currentprime
        allnums(notprime) = currentprime;
        count = count + 1;
        notprime = count*currentprime;
    end
    for i =currentprime+1:maxnumber+1  % Find the next prime
        if (allnums(i)==0)
            currentprime = i;
            break
        end
    end
end
nprimes = 0;
for i=2:maxnumber % Extract all the primes from allnums
    if (allnums(i)==0)
        nprimes = nprimes + 1;
        primes(nprimes) = i;
    end
end
end