A Shooting Method Approach to Boundary-Value ODEs

Topic: Differential Equations Boundary-Value Problems

Subject: Numerical Methods

Tool: Scilab

By: Gani Comia, July 2025

• The Shooting Method

The shooting method is a numerical technique of solving boundary-value problems (BVPs) for second-order ordinary differential equations (ODE) of the form shown in Equation (1).

$$y''(x) = f(x, y, y') , \quad a \leq x \leq b \tag{1}$$

with boundary conditions

$$y(a) = \alpha , \quad y(b) = \beta \tag{2}$$

A sample ODE, along with the boundary conditions provided in Equations (3) and (4), will be solved to demonstrate the method in this article.

$$\frac{d^2 y}{dx^2} + e^{-x y} + \sin\left(\frac{dy}{dx}\right) = 0 , \quad 1 \leq x \leq 2 \tag{3}$$

$$y(1) = 0 , \quad y(2) = 0 \tag{4}$$

The idea behind the shooting method is to convert the BVP to the initial-value problem (IVP) as in Equation (5) and (6). The given 2^nd-order ODE requires to be converted to the system of ODEs.

$$y'' = f(x, y, y') \Rightarrow \text{Sys of ODE} \left\{ \begin{aligned} y_1 &= y \\ y_2 &= y' \\ y_1^{'} &= y_2 \\ y_2^{'} &= f(x, y_1, y_2) \end{aligned} \right\} \tag{5}$$

$$y_1(a) = \alpha , \quad y_2(b) = s \tag{6}$$

Solving the converted IVP requires providing a guess for the value of \(s\) and solve the system numerically using Runge-Kutta method to get the value of the solution, \(y_{\beta} = y(b)\), at \(x = b\). If \(y_{\beta} \neq \beta\), adjust the guess, \(s\), so that \(y_{\beta} - \beta = 0\). For the given two guesses of \(s\), compute the best guess to make \(y_{\beta} - \beta = 0\) by employing the root-finding algorithm which is the Bisection method for this article.

For our given problem in Equation (3) and (4), two guesses value for \(s\) will be tried for \(y_2(a)\).

$$s_1 = 0.1 , \quad s_2 = \frac{\pi}{2} \tag{7}$$

Figure 1 is the resulting solution for the IVP using \(s_1\) and \(s_2\).

Figure 1. Solution to ODE IVP with \(y(a) = \alpha\) , \(s_1\) , and \(s_2\).

As shown in Figure 1, the two guesses, \(s_1\) and \(s_2\), do not meet the \(y(b) = \beta\) or \(y(2) = 0\) boundary conditions. Finding the value of \(s\) that makes \(y_{\beta} - \beta = 0\) requires defining a Residual function in Equation 8.

$$\phi(s) = y(b; s) - \beta \tag{8}$$

Residual function measures how far the solution with guesses is from satisfying the boundary condition at \(x = b\). For the given sample problem, \(\phi(s) = 1 \times 10^{-6}\) will be utilized in the root-finding algorithm, the Bisection method. Figure 2 is the visualization of the solution for sample ODE BVP.

Figure 2. Solution to ODE BVP using the Shooting Method.

The solution \(y(x)\) satisfies the boundary conditions \(y(a) = \alpha\) and \(y(b) = \beta\). Implementation of the Runge-Kutta method and the Bisection method is used in our Scilab script. The two guesses, \(s_1\) and \(s_2\), and the Residual function \(\phi(s)\) are used as parameters in the script.

• Scilab Script for Figure 2

// Copyright (C) 2025 - I.S.Comia
// Date of creation: Jul 12, 2025
// Shooting Method in Solving 2nd-Order ODE BVP
clear;clc;clf;
// Bundary Value Problem
// y'' = -exp(-x.*y)-sin(y'), y(1) = 0, y(2) = 0
// y'' = f(x,y,y'), y(a) = alpha, y(b) = beta

// Subroutine
// y'' = -exp(-x.*y)-sin(y') in Scilab function
function dydx=ode_system(x, y)
    // y(1) = y, y(2) = y'
    dydx = [y(2); -exp(x.*y(1))-sin(y(2))];
endfunction

// solving the ODE using IVP with guessed initial slope
// using Runge-Kutta Method
function yb=solve_ivp(s, a, b, h)
    // Initial value for y and y'
    x = a:h:b;
    n = length(x);
    y = zeros(2, n);
    y(:, 1) = [0; s];  // y(a) = 0, y'(a) = s (Note: s - slope)
    // this is Runge-Kutta
    for i = 1:n-1
        k1 = h * ode_system(x(i), y(:, i));
        k2 = h * ode_system(x(i) + h/2, y(:, i) + k1/2);
        k3 = h * ode_system(x(i) + h/2, y(:, i) + k2/2);
        k4 = h * ode_system(x(i) + h, y(:, i) + k3);
        y(:, i+1) = y(:, i) + (k1 + 2*k2 + 2*k3 + k4)/6;
    end

    yb = y(1, $); // return y(b)
endfunction

// Main Function
// Parameters
a = 1;
b = 2;
alpha = 0;
beta = 0;
h = 0.05;

// Two Initial guesses for s = y'(a) or slopes, s @ x = a
s1 = 0.1;
s2 = %pi/2;
tol = 1e-6;

// Plot of IVP using s1 and s2
x = a:h:b
x0 = a
y1_0 = [0;s1]
y1 = ode(y1_0,x0,x,ode_system)
y2_0 = [0;s2]
y2 = ode(y2_0,x0,x,ode_system)

plot(x, y1(1, :), "g.-", "linewidth",1);
plot(x, y2(1, :), "g.-", "linewidth",1);

// Shooting Method with Bisection Method
while abs(s2 - s1) > tol
    smid = (s1 + s2)/2;
    yb1 = solve_ivp(s1, a, b, h);
    ybmid = solve_ivp(smid, a, b, h);
    
    if (yb1 - beta)*(ybmid - beta) < 0 then
        s2 = smid;
    else
        s1 = smid;
    end
end

// Final solution with best guess
s_final = (s1 + s2)/2;
x = a:h:b;
n = length(x);
y = zeros(2, n);
y(:, 1) = [alpha; s_final];

for i = 1:n-1
    k1 = h * ode_system(x(i), y(:, i));
    k2 = h * ode_system(x(i) + h/2, y(:, i) + k1/2);
    k3 = h * ode_system(x(i) + h/2, y(:, i) + k2/2);
    k4 = h * ode_system(x(i) + h, y(:, i) + k3);
    y(:, i+1) = y(:, i) + (k1 + 2*k2 + 2*k3 + k4)/6;
end

// Visualization of solution to ODE BVP
fig = gcf()
fig.figure_size = [700, 700]
plot(x, y(1, :), "b-", "linewidth",3);
plot(x, y(1, :), "ro", "linewidth",4.5);
title("$\text{Shooting Method Solution to ODE BVP}$","fontsize",4.5);
xlabel("$\LARGE x$")
ylabel("$\LARGE y(x)$")
note = ["$\Large y(x) \; \text{w/ IVPs}$";"$\Large y(x) \; \text{w/ BVP}$";" "]
legend(note,2,with_box=%F)
xgrid(color("grey"),1,7)
// Plot annotation
h1 = xstring(1.55,0.725,"https://gani-mech-toolbox.blogspot.com")
h1.font_size = 2
eqn1 = "$\frac{d^2 y}{dx^2}+e^{-xy}+\sin\left(\frac{dy}{dx}\right)=0\text{,}\;

y(1)=0\text{,}\;y(2)=0$"
h2 = xstring(1.05,-0.4,eqn1)
h2.font_size = 3.5
h3 = xstring(1.05,-0.25,"Given : 2nd-Order ODE BVP")
h3.font_size = 3
eqn2 = "$y(a) = \alpha, \quad \left.\frac{dy}{dx} \right|_{x=a} = s_1$"
h4 = xstring(1.6,-0.125,eqn2)
h4.font_size = 3.5
eqn3 = "$y(a) = \alpha, \quad \left.\frac{dy}{dx} \right|_{x=a} = s_2$"
h5 = xstring(1.5,0.55,eqn3)
h5.font_size = 3.5
eqn4 = "$y(a) = \alpha, \quad \left.\frac{dy}{dx} \right|_{x=a} \; \mid \; 

y_\beta - \beta = 0$"
h6 = xstring(1.5,0.15,eqn4)
h6.font_size = 3.5

ax = gca()
ax.data_bounds = [1 -0.4; 2 0.8]

• Scilab Script for Figure 1

// Copyright (C) 2025 - I.S.Comia
// Date of creation: Jul 12, 2025
// Solution to 2nd-Order ODE IVP using Scilab's ode() function
clear;clc;clf;
// Initial Value Problem
// y'' = -exp(-x.*y)-sin(y'), y(1) = 0, y'(1) = s1 & s2
// y'' = f(x,y,y'), y(a) = alpha, y'(a) = s1 & s2

// Subroutine
// y'' = -exp(-x.*y)-sin(y') in Scilab function
function dydx=ode_system(x, y)
    // y(1) = y, y(2) = y'
    dydx = [y(2); -exp(x.*y(1))-sin(y(2))];
endfunction

// Main Function
// Parameters
a = 1;
b = 2;
alpha = 0;
h = 0.05;
// Two Initial guesses for s = y'(a) or slope s @ x = a
s1 = 0.1;
s2 = %pi/2;

// Plot of IVP using s1 and s2
x = a:h:b
x0 = a
y1_0 = [0;s1]
y1 = ode(y1_0,x0,x,ode_system)
y2_0 = [0;s2]
y2 = ode(y2_0,x0,x,ode_system)

// Visualization of solution to ODE IVP
fig = gcf()
fig.figure_size = [700, 700]
plot(x, y1(1, :), "g.-", "linewidth",1);
plot(x, y2(1, :), "m.-", "linewidth",1);
note = ["Solution to ODE IVP with ,","$y(a)=\alpha$",",","$\frac{dy}{dx}(a) = s$"]
title(note,"fontsize",3.5);
xlabel("$\LARGE x$")
ylabel("$\LARGE y(x)$")
note1 = ["$y(x) \; \text{w/ IVP} \; s_1 = 0.1$","$y(x) \; \text{w/ IVP} \; 

s_2 = \frac{\pi}{2}$"]
a1 = legend(note1,2,with_box=%F)
a1.font_size = 3
xgrid(color("grey"),1,7)
// Plot annotations
h1 = xstring(1.55,0.725,"https://gani-mech-toolbox.blogspot.com")
h1.font_size = 2
eqn1 = "$\frac{d^2 y}{dx^2}+e^{-xy}+\sin\left(\frac{dy}{dx}\right)=0\text{,}\;

y(1)=0\text{,}\;\frac{dy}{dx}(1)=s$"
h2 = xstring(1.05,-0.4,eqn1)
h2.font_size = 3.5
h3 = xstring(1.05,-0.25,"Given : 2nd-Order ODE IVP")
h3.font_size = 3
eqn2 = "$y(a) = \alpha, \quad \left.\frac{dy}{dx} \right|_{x=a} = s_1$"
h4 = xstring(1.6,-0.125,eqn2)
h4.font_size = 3.5
eqn3 = "$y(a) = \alpha, \quad \left.\frac{dy}{dx} \right|_{x=a} = s_2$"
h5 = xstring(1.5,0.55,eqn3)
h5.font_size = 3.5
// Plot limits
ax = gca()
ax.data_bounds = [1 -0.4; 2 0.8]

The shooting method is an effective approach for solving second-order boundary value problems (BVPs), particularly in mechanical, thermal, and structural engineering applications. To apply this method, the BVP is converted into an initial value problem (IVP) by estimating the necessary initial conditions. The IVP is then solved, and the results are compared with the boundary conditions at the opposite end. This process is repeated, adjusting the initial guesses as needed, until the solution satisfies the boundary conditions. While this approach involves some trial and error, linear second-order problems typically require no more than two iterations to converge.

Feel free to comment for inquiry, clarification, correction or suggestion for improvement. Drop your email to make a request to the author.

References

1. Gerald, C. and Wheatly, P., Applied Numerical Analysis. Addison-Wesley Publishing Company, Inc. 1994.
2. ChatGPT. AI Tools. Accessed on July 12, 2025.
3. Google AI Overview. AI Tools. Accessed on July 12, 2025.