function H2O2_kinetics()

time = [0;120;300;600;1200;1800;2400;3000;3600]; %time in seconds
data = [1;0.91;0.78;0.59;0.37;0.22;0.13;0.08;0.05]; %H2O2 conc in mol/L

%--------------------------------------------------------------------------
%FIRST-ORDER:
%--------------------------------------------------------------------------

%If this is a first-order reaction, this plot looks linear:
    log_of_data = log(data);

    figure(1)
    plot(time, log_of_data, 'b*','Markersize',8,'LineWidth',3)
    xlabel('time, in seconds')
    ylabel('log of concentration')
    title('Log of Concentration of H2O2')
    grid on

%Get an initial guess for the first-order rate constant k:
    coeffs1 = polyfit(time, log_of_data, 1);  %This gets us the best-fit
                                              %line.  The 1 indicates a
                                              %first-degree polynomial.
                                              %The variable named coeffs
                                              %will be a vector that
                                              %contains the slope and
                                              %y-intercept of the
                                              %best-fit line.
    slope_of_line = coeffs1(1);
    yint_of_line = coeffs1(2);
    
    initial_k_1st_order = -slope_of_line; %initial k is the negative of the 
                                          %slope of the best-fit line FOR
                                          %FIRST-ORDER model. See
                                          %Equation (2).
    
%Now look for the best parameters y0 and k:
    %We need a starting point for our search:
        init_params_1st_order = [data(1),initial_k_1st_order]; 
                             %= [y(0), k]
                                        
    %Now the big event!  The next command starts with our initial guess at
    %the parameters for the model, and returns the parameters that produce
    %the best fit for the data (best fit here means that it minimizes the
    %sum of squared errors).
        bestfit_params_1st_order ...
            = fminsearch(@(params) sum_of_squared_errors(params, ...
                                                      @first_order_model, ...
                                                      data, time), ...
                                                    init_params_1st_order);
                                                
    %bestfit_params_1st_order   %Uncomment this line if you want to
                                %print the parameters to the screen
                                
    %What does our best-fit model predict for the given times in the table?                            
        tvec = time(1) : 0.1 : time(end);
        bestfit_model_1st_order_output ...
                    = first_order_model(bestfit_params_1st_order, tvec);
        
    %Plot:
        figure(2)
        plot(time, data, 'b*', tvec, bestfit_model_1st_order_output,'k', ...
                'Markersize',8,'LineWidth',3)
        xlabel('time, in seconds')
        ylabel('concentration, in mol/L')
        title('Concentration of H2O2')
        legend('data','best-fit first-order model')
        grid on
%--------------------------------------------------------------------------
%SECOND-ORDER:
%--------------------------------------------------------------------------
%You'll probably want to begin here...

%--------------------------------------------------------------------------
%COMPARE SSEs (smallest SSE is an indicator of best model)
%--------------------------------------------------------------------------


end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Functions:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = sum_of_squared_errors(params, model, data, time)
    error = (model(params,time) - data);
    squared_errors = error.^2;
    out = sum(squared_errors);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = first_order_model(params,time)
    y0 = params(1);
    k = params(2);
    out = y0 * exp(-k * time);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%