clc, clear, close all;
data = readtable('ph_eq_data.xlsx', 'VariableNamingRule','preserve');
data = table2array(data);

set = [5e8; 5e9; 5e10; 5e11; 5e12; 5e13];
derivative = zeros(size(set,1),1);
integral = zeros(size(set,1),1);

for ci = 1:size(set,1)
    middle = ci * 2 + 1;
interpolatingdata = [data(middle - 2,:);data(middle - 1,:);data(middle,:);data(middle + 1,:)];

syms x;

A = [ones(4,1),interpolatingdata(:,1),interpolatingdata(:,1).^2,interpolatingdata(:,1).^3];
b = [interpolatingdata(:,2)];
X = linsolve(A,b);

syms x;

y = [1,x, x^2, x^3] * X;

yp = diff(y);

ym = matlabFunction(y);
ypm = matlabFunction(yp);
derivative(ci) = ypm(set(ci));
integral(ci) = ym(set(ci));
end

figure(1);
hold on;
plot(set, derivative)
