Execute MATLAB/Octave Online

z = [0 2.3 4.9 9.1 13.7 18.3 22.9 27.2]
T = [22.8 22.8 22.8 20.6 13.9 11.7 11.1 11.1]
z = z'
T = T'

A = bsxfun( @power, z, 0:3 )
b = A'*T
AtA = A' * A

fprintf("Doing LU....\n\n")
[m,n]=size(AtA);
U=zeros(m);
L=zeros(m);
for j=1:m
    L(j,j)=1;
end
for j=1:m
    U(1,j)=AtA(1,j);
end
for i=2:m
    for j=1:m
        for k=1:i-1
            s1=0;
            if k==1
                s1=0;
            else
            for p=1:k-1
                s1=s1+L(i,p)*U(p,k);
            end
            end
            L(i,k)=(AtA(i,k)-s1)/U(k,k);
           end
         for k=i:m
             s2=0;
           for p=1:i-1
               s2=s2+L(i,p)*U(p,k);
           end
           U(i,k)=AtA(i,k)-s2;
         end
    end
end
L
U
% [L,U] = lu(AtA);
% L([1 4],:) = L([4 1],:);
% L([2 3],:) = L([3 2],:)
% U

fprintf("Doing LDLt....\n\n")
D = [U(1, 1) 0 0 0
0 U(2, 2) 0 0
0 0 U(3, 3) 0
0 0 0 U(4, 4)]
L
Lt = L'

fprintf("Substitution.....\n\n")
y = L \ b
z = D \ y
x = Lt \ z

first_derivative_of_x = fliplr(x')
%in C/m
fprintf('The gradient is: %f + x(%f + x(%f + x(%f)))\n\n', first_derivative_of_x(4), first_derivative_of_x(3), first_derivative_of_x(2), first_derivative_of_x(1))

second_derivative_of_x = polyder(polyder(first_derivative_of_x))
inflection_points = roots(second_derivative_of_x)


flux = polyval(first_derivative_of_x, inflection_points)/100 * 0.01 * 86400 * -1

fprintf("\n\n================================== VANDERMONDE ==================================\n")

A = bsxfun( @power, z, 0:7 )
AtA = A' * A