A couple of notes. Since the integrand is never less than $1$, you know that $x_l<170$. In fact, when $x=170$, the integrand will be about $1.087$, so the eventual answer is bounded below by $\frac{170}{1.087}=156.5$. That gives you a pretty good idea about the starting point.
Bisection is useful when you don't have the derivative available, but here the you have the derivative in the form of the integrand, so it seems to be simpler to implement Newton's method in this case unless bisection is the point of the assignment.
Also you know the integral up to a point near to the actual root at each step after the first, so you need not integrate all the way from zero to the next guess for the root, just from the previous guess. Even if you integrated from zero every time, this program would zip through pretty fast anyway.
EDIT: Newton's method
% simp.m
f = @(x) sqrt(1+(x^2/68000)^2);
err = 1;
tol = 1.0e-8;
N = 10000;
a = 0;
b = 163;
while abs(err) > tol,
h = (b-a)/N;
y = f(a);
for i = 1:N/2,
y = y+4*f((2*i-1)*h)+2*f(2*i*h);
end
y = y+4*f(b-h)+f(b);
y = h/3*y-170;
yp = f(b);
err = y/yp;
b = b-err;
end
b
Once the the quadrature formula obtained by cubic interpolation of $f$ based at nodes $(a,f(a))$, $(a+h,f(a+h))$, $(a+2h,f(a+2+2h))$ and $(b,f(b))$ with $h=\frac{b-a}{3}$ and its error are established, notice that this formula can be rewritten as
$$\int^b_af\approx\frac{3h}{8}\Big(f(a)+3f(a+h)+3f(a+2h)+f(b)\Big)$$
with error $e_h=-\frac{(b-a)^5}{6480}f^{(iv)}(c)$ for some $c\in (a, b)$.
Divide the interval $[a,b]$ in $3m$ pieces of same length $h=\frac{b-a}{3m}$. This defines a partition $x_k=a+\frac{b-a}{3m}$.
On the each subinterval of the form $[x_{3k},x_{3(k+1)}]$, $k=0,\ldots ,m-1$ apply Simpson's rules based on cubic interpolation (the fist formula in the OP). That is
$$\begin{align}
\int^b_af=\sum^{m-1}_{k=0}\int^{x_{3(k+1)}}_{x_{3k}}f\tag{1}\label{one}
\end{align}$$
and
$$I_k:=\int^{x_{3(k+1)}}_{x_{3k}}f\approx\frac{3h}{8}\Big(f(x_{3k})+3f(x_{3k+1})+3f(x_{3k+2})+f(x_{3(k+1)})\Big)=A_k$$
Denoting the righthand side of \eqref{one} by $A$ we have that
\begin{align}
\frac{8}{3h}A&=\sum^{m-1}_{k=0}\Big(f(x_{3k})+3f(x_{3k+1})+3f(x_{3k+2})+f(x_{3(k+1)})\Big)\\
&=f(x_0)+\sum^{m-1}_{k=0}3\big(f(x_{3k+1})+f(x_{3k+2})\big) +2\sum^{m-1}_{k=1}f(x_{3k}) + f(x_{3m})
\end{align}
The composite formula follows.
As for the error, each approximation $A_k$ to $I_k$ has error $-\frac{3}{80} h^5 f^{(iv)}(c_k)$ for some $c_k\in (x_{3k},x_{3(k+1)})$. The sum of all errors is
$$E:=-\frac{3}{80}h^5\sum^{m-1}_{j=0}f^{(iv)}(c_k)$$
If we assume that $f\in\mathcal{C}^4([a,b])$, there there is $c\in(a, b)$ such that $f^{(iv)}(c)=\frac1m\sum^{m-1}_{k=0}f^{(iv)}(c_k)$ and so
$$
E=-\frac{1}{80}(b-a)h^4f^{(iv)}(c)
$$
With regards to the error in this quadrature, there are somewhat explicit formulas for the errors resulting from interpolation methods on evenly distributed nodes:
For $n\in\mathbb{N}$, $n\geq2$, the error of approximation to $\int^b_af$ in the quadrature method based on interpolating over the uniform partition $a_k=\frac{k}{n}(b-a)$, $0\leq k\leq n$, is
- If $n$ is even and $f\in C^{(n+2)}[a,b]$,
$$E_n=\frac{K_n}{(n+2)!}f^{(n+2)}(\eta),\quad K_n=\int^b_a x\prod^n_{k=0}(x-x_k)\,dx$$
for some $a<\eta<b$,
- If $n$ is odd and $f\in C^{(n+1)}[a,b]$,
$$E_n=\frac{K_n}{(n+1)!}f^{(n+1)}(\eta),\qquad K_n=\int^b_a\prod^n_{j=1}(x-x_j)\,dx$$
For Simpson's 3/8 rule, $n=3$ and so
$K_3=\int^b_a (x-a)(x-\tfrac{2a+
b}{3})(x-\tfrac{a+2b}{3})(x-b)\,dx$
The simplification, though messy with pencil and paper, yields
$E_3=-\frac{(b-a)^5}{6480}f^{(4)}(\eta)$ for some $a<\eta<b$.
Reference: Krylov, V. I. (translated by Stroud, A.), Approximate calculation of integrals, Dover Publications. 2005 pp 89-91
Best Answer
The formula of the Composite Simpson's 3/8 Rule is $$ \int_a^b{f(x) dx} \approx \frac{3h}{8} \left[ f(x_0) + 3 \sum_{i=1}^{m}{\left(f(x_{3i-2})+f(x_{3i-1})\right)} + \ 2 \sum_{i=1}^{m-1}f(x_{3i}) + f(x_{3m}) \right] $$ where $$ h=\frac{b-a}{3m}, $$ and for $i=0,1,2,\cdots,3m$ $$ x_i=a+ih. $$
Example.