npts = 20;
xlow = 2;
xhigh = 3.75;
x = linspace(xlow, xhigh, npts)';
fx = f(x);
x0 = 3;
fx0 = f(x0);
dfx0 = df(x0);
x1 = x0 - fx0/dfx0;
fx1 = f(x1);
gap = 0.75;
xzero = [xlow; xhigh];
yzero = [0; 0];
xt0 = [x1; x0 + gap];
yt0 = fx0 + dfx0 * (xt0-x0);
curvpts = [x0 fx0 x1 fx1];
axispts = [x0, 0, x1, 0;
x0, fx0, x1, fx1;
xlow, fx0, xlow, fx1];
tangents = [xzero, yzero, xt0, yt0];
plot(x, fx, axispts(:,1), axispts(:,2), axispts(:,3), axispts(:,4), ...
xzero, yzero, xt0, yt0)
table1 = [x fx];
save newton.dat table1 axispts tangents curvpts