function ethyl_alcohol_kinetics_kmb()

time = [0;393;699;1010;1265];
        %time, in seconds

data = [0; 0.0772; 0.1171; 0.1525; 0.1759];
        %concentration of ethyl alcohol, in mol/L

%--------------------------------------------------------------------------
%We want to find the k for the best-fit model.  We need a starting pnt:
    initial_k = 0.005;

%And we were asked to find the initial SSE in question 6a:
    initial_SSE = sum_of_squared_errors(initial_k, @model, data, time)

%Now we use fminsearch to find the k so that the model best fits the
%data:
    bestfit_k = fminsearch(@(k) sum_of_squared_errors(k,@model,data,time), ...
                                            initial_k)

%Let's find the SSE associated with the best-fit model (to confirm it's
%less than it was with our sub-optimal initial guess for k):
    bestfit_SSE = sum_of_squared_errors(bestfit_k, @model, data, time)

%We want to plot the output of our model with the data:
    tvec = time(1):0.1:time(end);
    bestfit_model_output = model(bestfit_k,tvec);

    figure(1)
    plot(time, data, 'b*', tvec, bestfit_model_output,'k', ...
            'Markersize',8,'LineWidth',3)
    xlabel('time, in seconds')
    ylabel('concentration, in mol/L')
    title('Concentration of the product Ethyl Alcohol')
    legend('data','best-fit model')
    grid on

end


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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = sum_of_squared_errors(k, model, data, time)
    error = (model(k,time) - data);
    squared_errors = error.^2;
    out = sum(squared_errors);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = model(k,time)
    A0 = 0.5638;
    B0 = 0.3114;
    out = (A0*B0*(exp(k*A0*time)-exp(k*B0*time)))...
                ./(A0*exp(k*A0*time)-B0*exp(k*B0*time));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%