function [t,C_a,C_b] = odesintegrate(solverfun,f0,tinterval,k1,k2);
%[t,C_a,C_b] = odesintegrate(solvername,f0,tinterval,k1,k2);
%GOALS:
%demonstrate the use of the various ode solvers to integrate a pair of
%interdependent odes
%
%show how to pass additional problem parameters to the function that
%calculates the derivatives 
%
%DESCRIPTION
%this example is a well stirred isothermal batch reactor with three species
%and two first order reactions:
%     k1  k2
%   A -> B -> C
%     
%only species A and B are followed as C can be gotten from a mass balance
%(C is not important and is not followed)
%
%VARIABLES
%solverfun: time integrator used, i.e. @ode45(explicit) or @ode15s(implicit)
%f0: initial values of function pair [C_a_0,C_b_0]
%tinterval: [0,tmax] integrates the functions from 0 to tmax
%           if more than 2 times are given, e.g. [0 t1 t2 t3 ...] then
%           only the results of specified timesteps are returned, however,
%           MATLAB may make extra "unseen" steps to maintain accuracy
%options: a list of options that set, e.g. the tolerance of the solver.
%         If [] is given, then the default options are used
%k1: scalar that affects rate of reaction 1
%k2: scalar that affects rate of reaction 2
%
%EXAMPLE sequence, ramping up k1 with constant k2
%[t,C_a,C_b] = odesintegrate(@ode45,[1,0],[0,6],1,1); length(t)
%[t,C_a,C_b] = odesintegrate(@ode45,[1,0],[0,6],10,1); length(t)
%[t,C_a,C_b] = odesintegrate(@ode45,[1,0],[0,6],100,1); length(t)
%[t,C_a,C_b] = odesintegrate(@ode45,[1,0],[0,6],1000,1); length(t)
%[t,C_a,C_b] = odesintegrate(@ode15s,[1,0],[0,6],1000,1); length(t)


%BEGIN CODE

%define param struct of k1 and k2
param.k1 = k1;
param.k2 = k2;

%Set tolerances (these are the default tolerances)
options = odeset('RelTol',10^-3,'AbsTol',10^-6);

%call solver, get back list of t values and f values
[t,fs] = solverfun(@odefun,tinterval,f0,options,param);
%specific calling sequences, without allowing for any solver
%%%%%%[t,fs] = ode45(@odefun,tinterval,f0,options,param);
%%%%%%[t,fs] = ode15s(@odefun,tinterval,f0,options,param);

%results include a column vector of t values and a matrix of f values
%the f matrix is really two horizontally concatenated f column vectors
%they can be separated using colon operations
C_a = fs(:,1); % <= C_a is the first column vector
C_b = fs(:,2); % <= C_b is the second column vector

%plot results
figure(1);

plot(t,C_a,t,C_b);
legend('C_a','C_b')
xlabel('t')
ylabel('f')
xlim([0,t(end)]);
ylim([0,max([max(C_a),max(C_b)])]);
title(['k1 = ',num2str(param.k1),' k2 = ',num2str(param.k2)])

return;

function derivs = odefun(t,fs,param);
%note that in the context of this function,
%a scalar t and a column vector [C_a;C_b] are passed to dfdt
%also, everything in param has been passed

%switch from list of f's to actual names for ease of formulation of
%derivative equations
C_a = fs(1);
C_b = fs(2);
k1 = param.k1;
k2 = param.k2;
%express derivatives in terms of C_a, C_b, and t
dC_adt = -k1*C_a;
dC_bdt =  k1*C_a - k2*C_b;

%put derivative results back in column vector format for MATLAB
derivs = [dC_adt;dC_bdt];
return;