% MATLAB file used in the solution of Activity 14.1 of 
% "A First Course on Kinetics and Reaction Engineering."

% Data provided in the problem statement, in consistent units
V = 750.; % mL
CA0 = 1.0; % mol/L
T = 70 + 273.15; % K
t = [1
2
3
4
5
6
7
8
9
10
12
15
20]; % min
CA = [0.874
0.837
0.800
0.750
0.572
0.626
0.404
0.458
0.339
0.431
0.249
0.172
0.185]; % mol/L

% calculate nA
nA = CA*V;

% fit 2nd order polynomial to the nA vs t data
p = polyfit(t,nA,2);
% use its derivative to estimate dnA/dt
dnAdt_poly = 2*p(1)*t + p(2);

% use forward differences to estimate dnA/dt for all but the last point
dnAdt_fd =zeros(length(nA),1);
for i = 1:length(nA)-1
    dnAdt_fd(i) = (nA(i+1)-nA(i))/(t(i+1)-t(i));
end

% use backward differences to estimate dnA/dt for all but the first point
dnAdt_bd =zeros(length(nA),1);
for i = 2:length(nA)
    dnAdt_bd(i) = (nA(i)-nA(i-1))/(t(i)-t(i-1));
end

% use central differences to estimate dnA/dt
dnAdt_cd =zeros(length(nA),1);
for i = 2:length(nA)-1
    dnAdt_cd(i) = (dnAdt_fd(i) + dnAdt_bd(i))/2;
end
dnAdt_cd(1) = dnAdt_fd(1);
dnAdt_cd(length(nA)) = dnAdt_bd(length(nA));

% fit using forward differences
disp('Forward Differences');
x = -nA(1:length(nA)-1);
y_hat = dnAdt_fd(1:length(nA)-1);
figure;
FitLinmSR;
title('Forward Differences');

% fit using backward differences
disp('Backward Differences');
x = -nA(2:length(nA));
y_hat = dnAdt_bd(2:length(nA));
figure;
FitLinmSR;
title('Backward Differences');

% fit using central differences
disp('Central Differences');
x = -nA;
y_hat = dnAdt_cd;
figure;
FitLinmSR;
title('Central Differences');

% fit using the second order polynomial
disp('Second Order Polynomial Approximation');
x = -nA;
y_hat = dnAdt_poly;
figure;
FitLinmSR;
title('Second Order Polynomial');