function HW1_solutions

% MAX SCORE 22 POINTS

close all  % This closes all figures

%% Problem 1

for i = 1:401
    x(i)=-10+(i-1)*20/400;
end

% [2 POINTS]

%% Problem 2
y = sin(5*x) + sin(5.5*x);  % You can do this with a loop too

figure1 = figure;
axes('Parent',figure1,'Position',[0.13 0.1286 0.7547 0.7571],'FontSize',14);
plot(x,y,'LineWidth',2,'Color',[1 0 0])
xlabel({'x'},'FontSize',18,'FontName','Arial');
ylabel({'y'},'FontSize',18,'FontName','Arial');
title({'Solution to problem 1'},'FontSize',24,'FontName','Arial');

% [2 POINTS]

% The fig labels, etc are not required

%% Problem 3

for i = 1:401
    v(i) = (i-1)*4*pi/400;
    if (sin(v(i))>cos(v(i)))
        w(i) = sin(v(i));
    else
        w(i) = cos(v(i));
    end
end

figure2 = figure;
axes('Parent',figure2,'Position',[0.13 0.1286 0.7547 0.7571],'FontSize',14);
plot(v,w,'LineWidth',2,'Color',[1 0 0])
xlabel({'v'},'FontSize',18,'FontName','Arial');
ylabel({'w'},'FontSize',18,'FontName','Arial');
title({'Solution to problem 2'},'FontSize',24,'FontName','Arial');

% [3 POINTS]


%% Problem 4

n1 = count_variables(x,0);
fprintf('No. variables exceeding 0 in x %4d \n', n1)

for i=1:17
   pivec(i) =  floor(pi*10^(i-1))-10*floor(pi*10^(i-2));
end

n2 = count_variables(pivec,7);
fprintf('No. of first 17 digits of pi exceeding 7 %4d \n',n2)

%[5 POINTS]

% Here's one solution to the challenge problem (not required)
pivec = compute_pi_digits(1000);
n2 = count_variables(pivec,7);
fprintf('No. of first 1000 digits of pi exceeding 7 %4d \n',n2)
% There are many other possible solutions, e.g. using MATLAB's variable
% precision arithmetic feature...

%% Solution to Problem 5
r=0.225;
K=1000000;
[times,N_vals] = ode45(@logistic_eq,[0,30],310384);

figure3 = figure;
axes('Parent',figure3,'Position',[0.13 0.1286 0.7547 0.7571],'FontSize',14);
plot(times,N_vals,'LineWidth',2,'Color',[1 0 0])
xlabel({'Time (years)'},'FontSize',18,'FontName','Arial');
ylabel({'Population(thousands)'},'FontSize',18,'FontName','Arial');
title({'US Population'},'FontSize',24,'FontName','Arial');

% [5 POINTS]


%% Solution to Problem 6
r=0.2;
mu=0.6;
[times,w_vals] = ode45(@bz_eq,[0,100],[1,1]);

figure5 = figure;
axes1 = axes('Parent',figure5,'Position',[0.13 0.1286 0.7547 0.7571],'FontSize',14);
plot1 = plot(times,w_vals,'LineWidth',2);
set(plot1(1),'Color',[1 0 0],'DisplayName','x');
set(plot1(2),'Color',[0 0 1],'DisplayName','y');
xlabel({'Time (sec)'},'FontSize',18,'FontName','Arial');
ylabel({'Concentration'},'FontSize',18,'FontName','Arial');
title({'Cubic autocatalator, r=0.2, \mu=0.6'},'FontSize',24,'FontName','Arial');
legend1 = legend(axes1,'show');
set(legend1,'Position',[0.7033 0.7472 0.1207 0.128]);

r=0.1;
mu=0.6;
[times,w_vals] = ode45(@bz_eq,[0,100],[1,1]);
figure6 = figure;
axes1 = axes('Parent',figure6,'Position',[0.13 0.1286 0.7547 0.7571],'FontSize',14);
plot1 = plot(times,w_vals,'LineWidth',2);
set(plot1(1),'Color',[1 0 0],'DisplayName','x');
set(plot1(2),'Color',[0 0 1],'DisplayName','y');
xlabel({'Time (sec)'},'FontSize',18,'FontName','Arial');
ylabel({'Concentration'},'FontSize',18,'FontName','Arial');
title({'Cubic autocatalator, r=0.1, \mu=0.6'},'FontSize',24,'FontName','Arial');
legend1 = legend(axes1,'show');
set(legend1,'Position',[0.7033 0.7472 0.1207 0.128]);

% [5 POINTS]


%% The ODE solver functions are below
    
    function dNdt = logistic_eq(t,N)
        dNdt = r*N*(K-N)/K;
    end

    function dwdt = bz_eq(t,w)
    x = w(1); y=w(2);
    dxdt = mu-x*y^2-r*x;
    dydt = x*y^2-y+r*x;
    dwdt = [dxdt;dydt];
    end


end
%% Assorted functions...
function number_of_vars = count_variables(vector,threshold)
% Brute force method for calculating the number of elements
% of a vector that exceed threshold

number_of_vars = 0;
for i = 1:length(vector)
    if (vector(i)>threshold)
        number_of_vars = number_of_vars + 1;
    end
end
    


end
function number_of_vars = count_variables_fancy(vector,threshold)
% Fancy method for calculating the number of elements
% of a vector that exceed threshold.  This doesnt even really need
% to go inside a function...

number_of_vars = sum((vector>threshold));
    
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    
