function EN40_HW1_Solutions

  %MAX TOTAL SCORE:  15 POINTS for this portion of hw
  %    Plus up to 8 points extra credit
  
  % Max score for mupad portion 18.
  % Total max score 33 + 8 points extra credit.
  
close all % This closes any open figure windows

%% Solution to Problem 8

vars = randi(10,[100,1]); % 100 random integers.  See help for details

%[2 POINTS]

number_of_fives = count_variables(vars,5)  % This counts the # of fives.
% Create 17 digits of pi (copied from HW assignment)
for i=1:17
   pivec(i) =  floor(pi*10^(i-1))-10*floor(pi*10^(i-2));
end
number_of_ones_in_pi = count_variables(pivec,1)


% [5 POINTS - 3 for the function, and 2 for using it correctly.
%  Solutions that use a built-in MATLAB function to do
%  the calculation should get the 2
%  points but not the credit for writing the function, which was
%  the point of the problem]
%

%% Solution to the optional_problem
%
% This tests the fancy counting function 
values = [1:1:10];  % This vector contains numbers from 1 to 10
number_of_numbers = count_variables_fancy(vars,values); % Count them
total_numbers = sum(number_of_numbers) % The total should be 100
%
% Now we test it on pi
%
no_digits = 100;
pivec = compute_pi_digits(no_digits);
values = [0:1:9];
pi_distribution = count_variables_fancy(pivec,values);

figure1 = figure;
axes1 = axes('Parent',figure1,'XTick',[0 1 2 3 4 5 6 7 8 9],'FontSize',14);
box(axes1,'on');
hold(axes1,'all');

% Create bar
bar(values,pi_distribution,'FaceColor',[1 0 0],...
    'EdgeColor',[0.0784313753247261 0.168627455830574 0.549019634723663],...
    'BarWidth',0.6);
xlabel({'Value of digit'},'FontSize',14);
ylabel({'Number of occurrences'},'FontSize',14);
tit = strcat('No. occurrences of digits in first  ',num2str(no_digits),...
' digits of \pi');
title(tit,'FontSize',16);


% [up to 5 POINTS Extra credit, depending on how sophisticated
%  the solution is]


%%  Problem 10

s = 0.05;
lambda = 0.2;
rvals = 0:0.005:0.2;
for i = 1:length(rvals)
  drdtplot(i) = economic_growth(0,rvals(i));
end

figure2 = figure;

axes1 = axes('Parent',figure2,'FontSize',14);
box(axes1,'on');
hold(axes1,'all');
plot(rvals,drdtplot);
xlabel({'Standard of Living'},'FontSize',14);
ylabel({'Growth rate'},'FontSize',14);
title('Growth rate in Solow model','FontSize',16);

[times1,sols1] = ode45(@economic_growth,[0,20],0.01);
[times2,sols2] = ode45(@economic_growth,[0,20],0.3);
figure3 = figure;
axes1 = axes('Parent',figure3,'FontSize',12);
box(axes1,'on');
hold(axes1,'all');
plot(times1,sols1,'LineWidth',2,'Color',[1 0 0]);
plot(times2,sols2,'LineWidth',2);
xlabel({'Time (years)'},'FontSize',14);
ylabel({'Living Standard'},'FontSize',14);
title({'Living standard predicted by Solow growth model'},'FontSize',16);

    function drdt = economic_growth(t,r)
       drdt = s*r^(1/3)-lambda*r;        
    end


% [2 points for the function economic_growth function
%  2 points for each plot]
% Nicely labeled plots, etc can get 1 point extra credit each

%% Problem 11

   N0 = 0.25;
   B0 = 0.1;
   beta = 0.5;
   alpha = 2.7;
   
   [times1,sols1] = ode45(@chemostat,[0,20],[B0;N0]);
   beta = 0.8;
   alpha = 2.7;
   [times2,sols2] = ode45(@chemostat,[0,20],[B0;N0]);
   figure4 = figure;
   axes1 = axes('Parent',figure4,'FontSize',14);
   box(axes1,'on');
   hold(axes1,'all');
   plot1 = plot(times1,sols1,'Parent',axes1,'LineWidth',2,'LineStyle','--');
   set(plot1(1),'Color',[1 0 0],'DisplayName','B, \beta=0.5');
   set(plot1(2),'Color',[0 0 1],'DisplayName','N, \beta=0.5');
   plot2 = plot(times2,sols2,'Parent',axes1,'LineWidth',2);
   set(plot2(1),'Color',[1 0 0],'DisplayName','B, \beta=0.8');
   set(plot2(2),'Color',[0 0 1],'DisplayName','N, \beta=0.8');
   xlabel('Time','FontSize',16);
   ylabel('Concentration','FontSize',16);
   legend1 = legend(axes1,'show');
   set(legend1,...
    'Position',[0.759997203579418 0.659403669724771 0.108221476510067 0.237003058103976]);
   
    function dwdt = chemostat(t,w)
        B = w(1); N = w(2); % Extract the variables from the vector
        dBdt = alpha*B*N/(1+N)-B; 
        dNdt = -B*N/(1+N) - N + beta;
        dwdt = [dBdt;dNdt];
    end

 % [2 points for the chemostat function
 %  2 points for each plot
 %  Nicely labeled plots can get 1 point extra credit each.

 
end
function number_of_vars = count_variables(vector,value)
     
  number_of_vars = 0;
  for i=1:length(vector)
      if (vector(i)==value)
          number_of_vars = number_of_vars+1;
      end
  end
     
        
end
function number_of_vars = count_variables_fancy(vector,values)
% This is a fancier function that allows a list of values
% to be provided - each one in the list will be counted.
  number_of_vars = zeros(length(values),1);
  
  for i=1:length(vector)
    for k=1:length(values)
      if (vector(i)==values(k))
          number_of_vars(k) = number_of_vars(k)+1;
      end
    end
  end
     
        
end
function pivec = compute_pi_digits(n)
% Here's the challenge problem - there are many clever ways to compute
% the digits of pi - this is an implementation of the Rabinowitz-Wagon
% spigot algorithm. 
% www.mathpropress.com/stan/bibliography/spigot.pdf
      if (n==1)
          pivec(1) = 3;
      else
        A=2*ones(1,ceil(10*n/3));
        r=0;
        for j=1:n
          A=A*10;
          for k=floor(10*n/3):-1:2
              A(k)=A(k)+r;
              r=floor(A(k)/(2*k-1))*(k-1);
              A(k)=mod(A(k),2*k-1);
          end
          A(1)=A(1)+r;
          value=floor(A(1)/10);
          A(1)=A(1)-10*value;
          pivec(j)=value;
        end
      end
end    