The main point of Shooting methods is to take a boundary value problem, which is typically difficult to solve, and replace it with an initial value problem, which is easier to solve.  In our case, we will solve the boundary value problem
\begin{eqnarray*}
y'' = y && y(-1) = y(1) = 1
\end{eqnarray*}
Which we then turn into the initial value problem
\begin{eqnarray*}
y'' = y & y(-1) = 1 & y'(-1) = \eta
\end{eqnarray*}
Where $\eta$ is unknown.  We have an implementation (in Python) that will compute the solution of this IVP using Euler's method (with $h = 1/50$), as a function of $\eta$.  In order to find a good value for $\eta$, we define the function
\[F(\eta) = y(1,\eta) - 1\]
This equation is the difference between the solution we want and the solution we get using Euler's method.  Then we want to solve for $0$ using bisection method.  (Remembering that $y(1,\eta)$ is found using our implementation of Euler's method.)  For bisection, we need to find two values that bracket the root of our function $F(\eta)$.  We do this, more or less, by hand.  To begin with, we iterate over the interval $[-1,1]$ by steps of $0.1$ and get the values
\begin{center}
\begin{tabular}{|l|l|} \hline
$\eta$ & $F(\eta)$ \\ \hline
-1.0 & -0.867380444105 \\ \hline
-0.9 & -0.511779115987 \\ \hline
-0.8 & -0.156177787869 \\ \hline
-0.7 & 0.199423540248 \\ \hline
-0.6 & 0.555024868366 \\ \hline
-0.5 & 0.910626196484 \\ \hline
-0.4 & 1.2662275246 \\ \hline
-0.3 & 1.62182885272 \\ \hline
-0.2 & 1.97743018084 \\ \hline
-0.1 & 2.33303150896 \\ \hline
0 & 2.68863283707 \\ \hline
0.1 & 3.04423416519 \\ \hline
0.2 & 3.39983549331 \\ \hline
0.3 & 3.75543682143 \\ \hline
0.4 & 4.11103814955 \\ \hline
0.5 & 4.46663947766 \\ \hline
0.6 & 4.82224080578 \\ \hline
0.7 & 5.1778421339 \\ \hline
0.8 & 5.53344346202 \\ \hline
0.9 & 5.88904479013 \\ \hline
\end{tabular}
\end{center}
To illustrate how shooting works, we take a set of solutions, say the ones for $\eta = -1,-0.9,\ldots,-0.5$ and plot them along with the numerical solution.  
\begin{center}
\includegraphics[height=9cm,width=9cm]{all_shooting}
\end{center}
We can use the values in our table to determine a bracket for bisection.  Using SciPy's bisection method on the bracket $-0.8,-0.7$ we, find the approximate value of $\eta = -0.756080651134$.  If we plot the solution given by Euler's method by itself, using $h = 1/50$, we get this plot
\begin{center}
\includegraphics[height=9cm,width=9cm]{numsol_shooting}
\end{center}
The true solution is given by
\[y(x) = \frac{e^{x+1} + e^{1-x}}{e^{2} + 1}\]
If we plot the true solution along with our approximated solution, we get
\begin{center}
\includegraphics[height=9cm,width=9cm]{true_numsol_shooting}
\end{center}
There is virtually no difference between the true solution and the approximated one.  It is interesting to note that we are essentially relying on the accuracy of both Euler's method and the bisection method, therefore, we don't have to do any extra error analysis as our results are as good as the methods we choose.  For example, we could use more accurate methods or different methods altogether.  (Code available upon request.)
