Add this code to the end of your script:
% Fit the Bex points
coeffs1 = polyfit(X, Bex, 2);
newX = linspace(X(1), X(end), 100);
fittedValues = polyval(coeffs1, newX);
hold on;
plot(newX, fittedValues, 'b-');
% Fit the Bth points
coeffs2 = polyfit(X, Bth, 2);
newX2 = linspace(X(1), X(end), 100);
fittedValues2 = polyval(coeffs2, newX2);
hold on;
plot(newX2, fittedValues2, 'r-');
This will fit a quadratic to your data points. You can use another order if you want. Do you know what model you want? If you don't want a polynomial, you'll have to use a different function. Look up fitting in the help. There is also a Curve Fitting Toolbox that you may have that could help.
