function HW1_solutions
%
% EN0040 Dynamics and Vibrations
% Solutions to Homework 1, Spring 2015
% Total: 15 POINTS
%
    close all
%
%%   Problem 8 [3 POINTS]
%
    x = [0:1/101:1];
    y = zeros(101,1);
    for i =1:length(x)
        y(i) = logistic_map(0.5,x(i));
    end
    figure1 = figure;
    axes1 = axes('Parent',figure1,'FontSize',14);
    box(axes1,'on');
    hold(axes1,'all');
    plot(x,y,'Parent',axes1,'LineWidth',2,'Color',[1 0 0])
    xlabel('x','FontSize',14);
    ylabel('y','FontSize',14);
 
%
    %% Problem 9 [4 POINTS]
%
    x = zeros(31,1);
    x(1) = 0.5;
    for i=1:30
        x(i+1) = logistic_map(1.75,x(i));
    end
    figure2 = figure;
    axes2 = axes('Parent',figure2,'FontSize',14);
    plot1 = plot(x,'Parent',axes2,'LineWidth',2,'Color',[1 0 0],'DisplayName','r = 1.75');

    x = zeros(31,1);
    x(1) = 0.5;
    for i=1:30
        x(i+1) = logistic_map(3.6,x(i));
    end
    hold on
    plot2 = plot(x,'Parent',axes2,'LineWidth',2,'Color',[0 0 1],'DisplayName','r = 3.6');
    xlim([0 30]);
    xlabel('n','FontSize',14);
    ylabel('x_n','FontSize',14);
    title('Iterated Logistic Map');
    legend1 = legend(axes2,'show');
    set(legend1,...
    'Position',[0.750880902043693 0.834786821705426 0.130021141649049 0.0784883720930233]);
%
%%  Problem 10 [4 POINTS]
%
    s = 0.817;   
    r = 731000;  % kg/year
    m = 116000;  % kg/year
    q = 7.88;    % Exponent (dimensionless)
%    
    initial_p = 75000;
    L = 38000;   % kg/year
%    
    tstart = 0;
    tstop = 10;
%    
    [time,p] = ode45(@diff_eq_prob_10,[tstart,tstop],initial_p);
    figure3 = figure;
    axes3 = axes('Parent',figure3,'FontSize',14);
    plot(time,p/1000,'LineWidth',2,'Color',[0 0 1],'DisplayName','L = 38000kg/s');
    hold on
%
    L=38500;   % kg/year
    [time,p] = ode45(@diff_eq_prob_10,[tstart,tstop],initial_p);
    plot(time,p/1000,'LineWidth',2,'Color',[1 0 0],'DisplayName','L = 38500kg/s');
    xlabel('time (years)','FontSize',14);
    ylabel('P (metric tonnes)','FontSize',14);
    title('Phosphorus Concentration');
    legend1 = legend(axes3,'show');
    set(legend1,...
    'Position',[0.274665257223398 0.785368217054263 0.190274841437632 0.0784883720930233]);
%
% Function for problem 10    
%
    function dpdt = diff_eq_prob_10(t,p)
       dpdt = L - s*p + r*(p^q/(m^q+p^q)); 
    end

%
%%  Problem 11 [4 POINTS]
%
    a = 0.1;
    b = 0.1;
    c = 18;
    
    tstart = 0;
    tstop = 300;
    
    initial_w = [1;1;1];
    
    options=odeset('RelTol',0.00001);  %This increases accuracy - it is optional
    [time,solution] = ode45(@diff_eq_prob_11,[tstart,tstop],initial_w);
    figure4 = figure;
    axes4 = axes('Parent',figure4,'FontSize',14);
    plot3(solution(:,1),solution(:,2),solution(:,3))
    xlabel('x','FontSize',14);
    ylabel('y','FontSize',14);    
    zlabel('z','FontSize',14);  
    title('Trajecotry for the Rossler Attractor');
    
    function dwdt = diff_eq_prob_11(t,w)
        x = w(1); y=w(2); z = w(3);
        dxdt = -y-z;
        dydt = x+a*y;
        dzdt = b + z*(x-c);
        dwdt = [dxdt;dydt;dzdt];
    end
       

end
%% Function for problems 8 and 9
function y = logistic_map(r,x)
   y = r*x*(1-x);
end
