Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_ode_bvp_shoot_genpar (d02hb)

Purpose

nag_ode_bvp_shoot_genpar (d02hb) solves a two-point boundary value problem for a system of ordinary differential equations, using initial value techniques and Newton iteration; it generalizes function nag_ode_bvp_shoot_bval (d02ha) to include the case where parameters other than boundary values are to be determined.

Syntax

[p, soln, w, ifail] = d02hb(p, pe, e, m1, fcn, bc, range, 'n1', n1, 'n', n)
[p, soln, w, ifail] = nag_ode_bvp_shoot_genpar(p, pe, e, m1, fcn, bc, range, 'n1', n1, 'n', n)

Description

nag_ode_bvp_shoot_genpar (d02hb) solves a two-point boundary value problem by determining the unknown parameters p1,p2,,pn1${p}_{1},{p}_{2},\dots ,{p}_{{\mathit{n}}_{1}}$ of the problem. These parameters may be, but need not be, boundary values; they may include eigenvalue parameters in the coefficients of the differential equations, length of the range of integration, etc. The notation and methods used are similar to those of nag_ode_bvp_shoot_bval (d02ha) and you are advised to study this first. (The parameters p1,p2,,pn1${p}_{1},{p}_{2},\dots ,{p}_{{\mathit{n}}_{1}}$ correspond precisely to the unknown boundary conditions in nag_ode_bvp_shoot_bval (d02ha).) It is assumed that we have a system of n$\mathit{n}$ first-order ordinary differential equations of the form:
 (dyi)/(dx) = fi(x,y1,y2, … ,yn),  i = 1,2, … ,n, $dyi dx =fi(x,y1,y2,…,yn), i=1,2,…,n,$
and that the derivatives fi${f}_{i}$ are evaluated by fcn. The system, including the boundary conditions given by bc and the range of integration given by range, involves the n1${\mathit{n}}_{1}$ unknown parameters p1,p2,,pn1${p}_{1},{p}_{2},\dots ,{p}_{{\mathit{n}}_{1}}$ which are to be determined, and for which initial estimates must be supplied. The number of unknown parameters n1${\mathit{n}}_{1}$ must not exceed the number of equations n$\mathit{n}$. If n1 < n${\mathit{n}}_{1}<\mathit{n}$, we assume that (nn1)$\left(\mathit{n}-{\mathit{n}}_{1}\right)$ equations of the system are not involved in the matching process. These are usually referred to as ‘driving equations’; they are independent of the parameters and of the solutions of the other n1${\mathit{n}}_{1}$ equations. In numbering the equations for fcn, the driving equations must be put first.
The estimated values of the parameters are corrected by a form of Newton iteration. The Newton correction on each iteration is calculated using a Jacobian matrix whose (i,j)$\left(i,j\right)$th element depends on the derivative of the i$i$th component of the solution, yi${y}_{i}$, with respect to the j$j$th parameter, pj${p}_{j}$. This matrix is calculated by a simple numerical differentiation technique which requires n1${\mathit{n}}_{1}$ evaluations of the differential system.
If the parameter ifail is set appropriately, the function automatically prints messages to inform you of the flow of the calculation. These messages are discussed in detail in Section [Further Comments].
nag_ode_bvp_shoot_genpar (d02hb) is a simplified version of nag_ode_bvp_shoot_genpar_algeq (d02sa) which is described in detail in Gladwell (1979).

References

Gladwell I (1979) The development of the boundary value codes in the ordinary differential equations chapter of the NAG Library Codes for Boundary Value Problems in Ordinary Differential Equations. Lecture Notes in Computer Science (eds B Childs, M Scott, J W Daniel, E Denman and P Nelson) 76 Springer–Verlag

Parameters

You are strongly recommended to read Sections [Description] and [Further Comments] in conjunction with this section.

Compulsory Input Parameters

1:     p(n1) – double array
n1, the dimension of the array, must satisfy the constraint 1n1n$1\le {\mathbf{n1}}\le {\mathbf{n}}$.
An estimate for the i$\mathit{i}$th parameter, pi${p}_{\mathit{i}}$, for i = 1,2,,n1$\mathit{i}=1,2,\dots ,{\mathit{n}}_{1}$.
2:     pe(n1) – double array
n1, the dimension of the array, must satisfy the constraint 1n1n$1\le {\mathbf{n1}}\le {\mathbf{n}}$.
The elements of pe must be given small positive values. The element pe(i)${\mathbf{pe}}\left(i\right)$ is used
 (i) in the convergence test on the i$i$th parameter in the Newton iteration, and (ii) in perturbing the i$i$th parameter when approximating the derivatives of the components of the solution with respect to this parameter for use in the Newton iteration.
The elements pe(i)${\mathbf{pe}}\left(i\right)$ should not be chosen too small. They should usually be several orders of magnitude larger than machine precision.
Constraint: pe(i) > 0.0${\mathbf{pe}}\left(\mathit{i}\right)>0.0$, for i = 1,2,,n1$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$.
3:     e(n) – double array
n, the dimension of the array, must satisfy the constraint n2${\mathbf{n}}\ge 2$.
The elements of e must be given positive values. The element e(i)${\mathbf{e}}\left(i\right)$ is used in the bound on the local error in the i$i$th component of the solution yi${y}_{i}$ during integration.
The elements e(i)${\mathbf{e}}\left(i\right)$ should not be chosen too small. They should usually be several orders of magnitude larger than machine precision.
Constraint: e(i) > 0.0${\mathbf{e}}\left(\mathit{i}\right)>0.0$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,{\mathbf{n}}$.
4:     m1 – int64int32nag_int scalar
A value which controls exit values.
m1 = 1${\mathbf{m1}}=1$
The final solution is not calculated.
m1 > 1${\mathbf{m1}}>1$
The final values of the solution at interval (length of range)/(m11)$\left({\mathbf{m1}}-1\right)$ are calculated and stored sequentially in the array soln starting with the values of the solutions evaluated at the first end point (see range) stored in the first column of soln.
Constraint: m11${\mathbf{m1}}\ge 1$.
5:     fcn – function handle or string containing name of m-file
fcn must evaluate the functions fi${f}_{\mathit{i}}$ (i.e., the derivatives yi${y}_{\mathit{i}}^{\prime }$), for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$, at a general point x$x$.
[f] = fcn(x, y, p)

Input Parameters

1:     x – double scalar
x$x$, the value of the argument.
2:     y( : $:$) – double array
yi${y}_{\mathit{i}}$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$, the value of the argument.
3:     p( : $:$) – double array
The current estimate of the parameter pi${p}_{\mathit{i}}$, for i = 1,2,,n1$\mathit{i}=1,2,\dots ,{\mathit{n}}_{1}$.

Output Parameters

1:     f( : $:$) – double array
The value of fi${f}_{\mathit{i}}$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$. The fi${f}_{i}$ may depend upon the parameters pj${p}_{\mathit{j}}$, for j = 1,2,,n1$\mathit{j}=1,2,\dots ,{\mathit{n}}_{1}$. If there are any driving equations (see Section [Description]) then these must be numbered first in the ordering of the components of f in fcn.
6:     bc – function handle or string containing name of m-file
bc must place in g1 and g2 the boundary conditions at a$a$ and b$b$ respectively (see range).
[g1, g2] = bc(p)

Input Parameters

1:     p( : $:$) – double array
An estimate of the parameter pi${p}_{\mathit{i}}$, for i = 1,2,,n1$\mathit{i}=1,2,\dots ,{\mathit{n}}_{1}$.

Output Parameters

1:     g1( : $:$) – double array
The value of yi(a)${y}_{\mathit{i}}\left(a\right)$, (where this may be a known value or a function of the parameters pj${p}_{\mathit{j}}$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$ and j = 1,2,,n1$\mathit{j}=1,2,\dots ,{\mathit{n}}_{1}$).
2:     g2( : $:$) – double array
The value of yi(b)${y}_{\mathit{i}}\left(b\right)$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$, (where these may be known values or functions of the parameters pj${p}_{\mathit{j}}$, for j = 1,2,,n1$\mathit{j}=1,2,\dots ,{\mathit{n}}_{1}$). If n > n1$\mathit{n}>{\mathit{n}}_{1}$, so that there are some driving equations, then the first nn1$\mathit{n}-{\mathit{n}}_{1}$ values of g2 need not be set since they are never used.
7:     range – function handle or string containing name of m-file
range must evaluate the boundary points a$a$ and b$b$, each of which may depend on the parameters p1,p2,,pn1${p}_{1},{p}_{2},\dots ,{p}_{{\mathit{n}}_{1}}$. The integrations in the shooting method are always from a$a$ to b$b$.
[a, b] = range(p)

Input Parameters

1:     p( : $:$) – double array
The current estimate of the i$\mathit{i}$th parameter, pi${p}_{\mathit{i}}$, for i = 1,2,,n1$\mathit{i}=1,2,\dots ,{\mathit{n}}_{1}$.

Output Parameters

1:     a – double scalar
a$a$, one of the boundary points.
2:     b – double scalar
The second boundary point, b$b$. Note that b > a${\mathbf{b}}>{\mathbf{a}}$ forces the direction of integration to be that of increasing x$x$. If a and b are interchanged the direction of integration is reversed.

Optional Input Parameters

1:     n1 – int64int32nag_int scalar
Default: The dimension of the arrays p, pe. (An error is raised if these dimensions are not equal.)
n1${\mathit{n}}_{1}$, the number of parameters.
Constraint: 1n1n$1\le {\mathbf{n1}}\le {\mathbf{n}}$.
2:     n – int64int32nag_int scalar
Default: The dimension of the array e.
n$\mathit{n}$, the total number of differential equations.
Constraint: n2${\mathbf{n}}\ge 2$.

sdw

Output Parameters

1:     p(n1) – double array
The corrected value for the i$i$th parameter, unless an error has occurred, when it contains the last calculated value of the parameter.
2:     soln(n,m1) – double array
The solution when m1 > 1${\mathbf{m1}}>1$.
3:     w(n,sdw) – double array
sdw3n + 14 + max (11,n)$\mathit{sdw}\ge 3{\mathbf{n}}+14+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(11,{\mathbf{n}}\right)$.
With ${\mathbf{ifail}}={\mathbf{2}}$, 3${\mathbf{3}}$, 4${\mathbf{4}}$ or 5${\mathbf{5}}$ (see Section [Error Indicators and Warnings]), w(i,1)${\mathbf{w}}\left(\mathit{i},1\right)$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,\mathit{n}$, contains the solution at the point x$x$ when the error occurred. w(1,2)${\mathbf{w}}\left(1,2\right)$ contains x$x$.
4:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
One or more of the parameters n, n1, m1, sdw, e or pe is incorrectly set.
ifail = 2${\mathbf{ifail}}=2$
The step length for the integration became too short whilst calculating the residual (see Section [Further Comments]).
ifail = 3${\mathbf{ifail}}=3$
No initial step length could be chosen for the integration whilst calculating the residual.
Note: ${\mathbf{ifail}}={\mathbf{2}}$ or 3${\mathbf{3}}$ can occur due to choosing too small a value for e or due to choosing the wrong direction of integration. Try varying e and interchanging a$a$ and b$b$. These error exits can also occur for very poor initial choices of the parameters in the array p and, in extreme cases, because nag_ode_bvp_shoot_genpar (d02hb) cannot be used to solve the problem posed.
ifail = 4${\mathbf{ifail}}=4$
As for ${\mathbf{ifail}}={\mathbf{2}}$ but the error occurred when calculating the Jacobian.
ifail = 5${\mathbf{ifail}}=5$
As for ${\mathbf{ifail}}={\mathbf{3}}$ but the error occurred when calculating the Jacobian.
ifail = 6${\mathbf{ifail}}=6$
The calculated Jacobian has an insignificant column. This can occur because a parameter pi${p}_{i}$ is incorrectly entered when posing the problem.
Note: ${\mathbf{ifail}}={\mathbf{4}}$, 5${\mathbf{5}}$ or 6${\mathbf{6}}$ usually indicate a badly scaled problem. You may vary the size of pe. Otherwise the use of the more general nag_ode_bvp_shoot_genpar_algeq (d02sa) which affords more control over the calculations is advised.
ifail = 7${\mathbf{ifail}}=7$
The linear algebra function used (nag_lapack_dgesvd (f08kb)) has failed. This error exit should not occur and can be avoided by changing the initial estimates pi${p}_{i}$.
ifail = 8${\mathbf{ifail}}=8$
The Newton iteration has failed to converge. This can indicate a poor initial choice of parameters pi${p}_{i}$ or a very difficult problem. Consider varying the elements pe(i)${\mathbf{pe}}\left(i\right)$ if the residuals are small in the monitoring output. If the residuals are large, try varying the initial parameters pi${p}_{i}$.
ifail = 9${\mathbf{ifail}}=9$
ifail = 10${\mathbf{ifail}}=10$
ifail = 11${\mathbf{ifail}}=11$
ifail = 12${\mathbf{ifail}}=12$
ifail = 13${\mathbf{ifail}}=13$
Indicates that a serious error has occurred. Check all array subscripts and function parameter lists in the call to nag_ode_bvp_shoot_genpar (d02hb). Seek expert help.

Accuracy

If the process converges, the accuracy to which the unknown parameters are determined is usually close to that specified by you; the solution, if requested, may be determined to a required accuracy by varying e.

The time taken by nag_ode_bvp_shoot_genpar (d02hb) depends on the complexity of the system, and on the number of iterations required. In practice, integration of the differential equations is by far the most costly process involved.
Wherever they occur in the function, the error parameters contained in the arrays e and pe are used in ‘mixed’ form; that is e(i)${\mathbf{e}}\left(i\right)$ always occurs in expressions of the form
 e(i) × (1 + |yi|) $ei×(1+|yi|)$
and pe(i)${\mathbf{pe}}\left(i\right)$ always occurs in expressions of the form
 pe(i) × (1 + |pi|). $pei×(1+|pi|).$
Though not ideal for every application, it is expected that this mixture of absolute and relative error testing will be adequate for most purposes.
You may determine a suitable direction of integration a$a$ to b$b$ and suitable values for e(i)${\mathbf{e}}\left(i\right)$ by integrations with nag_ode_ivp_rkts_range (d02pe). The best direction of integration is usually the direction of decreasing solutions. You are strongly recommended to set ifail to obtain self-explanatory error messages, and also monitoring information about the course of the computation. You may select the channel numbers on which this output is to appear by calls of nag_file_set_unit_error (x04aa) (for error messages) or nag_file_set_unit_advisory (x04ab) (for monitoring information) – see Section [Example] for an example. Otherwise the default channel numbers will be used. The monitoring information produced at each iteration includes the current parameter values, the residuals and 2$2$-norms: a basic norm and a current norm. At each iteration the aim is to find parameter values which make the current norm less than the basic norm. Both these norms should tend to zero as should the residuals. (They would all be zero if the exact parameters were used as input.) For more details, in particular about the other monitoring information printed, you are advised to consult the specification of nag_ode_bvp_shoot_genpar_algeq (d02sa), and especially the description of the parameter monit there.
The computing time for integrating the differential equations can sometimes depend critically on the quality of the initial estimates for the parameters pi${p}_{i}$. If it seems that too much computing time is required and, in particular, if the values of the residuals printed by the monitoring function are much larger than the expected values of the solution at b$b$, then the coding of fcn, bc and range should be checked for errors. If no errors can be found, an independent attempt should be made to improve the initial estimates for pi${p}_{i}$.
The function can be used to solve a very wide range of problems, for example:
 (a) eigenvalue problems, including problems where the eigenvalue occurs in the boundary conditions; (b) problems where the differential equations depend on some parameters which are to be determined so as to satisfy certain boundary conditions (see Example 2 in Section [Example]); (c) problems where one of the end points of the range of integration is to be determined as the point where a variable yi${y}_{i}$ takes a particular value (see Example 2 in Section [Example]); (d) singular problems and problems on infinite ranges of integration where the values of the solution at a$a$ or b$b$ or both are determined by a power series or an asymptotic expansion (or a more complicated expression) and where some of the coefficients in the expression are to be determined (see Example 1 in Section [Example]); and (e) differential equations with certain terms defined by other independent (driving) differential equations.

Example

For this function two examples are presented. There is a single example program for nag_ode_bvp_shoot_genpar (d02hb), with a main program and the code to solve the two example problems given in Example 1 (EX1) and Example 2 (EX2).
Example 1 (EX1)
This example finds the solution of the differential equation
 y ′ ′ = (y3 − y′) / 2x $y′′=(y3-y′)/2x$
on the range 0x16$0\le x\le 16$, with boundary conditions y(0) = 0.1$y\left(0\right)=0.1$ and y(16) = 1 / 6$y\left(16\right)=1/6$. We cannot use the differential equation at x = 0$x=0$ because it is singular, so we take a truncated power series expansion
 y(x) = 1 / 10 + p1 × sqrt(x) / 10 + x / 100 $y(x)=1/10+p1×x/10+x/100$
near the origin where p1${p}_{1}$ is one of the parameters to be determined. We choose the interval as [0.1,16]$\left[0.1,16\right]$ and setting p2 = y(16)${p}_{2}={y}^{\prime }\left(16\right)$, we can determine all the boundary conditions. We take X1 = 16$\mathrm{X1}=16$. We write y = y(1)$y={\mathbf{y}}\left(1\right)$, y = y(2)${y}^{\prime }={\mathbf{y}}\left(2\right)$, and estimate PARAM(1) = 0.2$\mathrm{PARAM}\left(1\right)=0.2$, PARAM(2) = 0.0$\mathrm{PARAM}\left(2\right)=0.0$. Note the call to nag_file_set_unit_advisory (x04ab) before the call to nag_ode_bvp_shoot_genpar (d02hb).
Example 2 (EX2)
This example finds the gravitational constant p1${p}_{1}$ and the range p2${p}_{2}$ over which a projectile must be fired to hit the target with a given velocity.
The differential equations are
 y′ = tanφ v′ = ( − (p1sinφ + 0.00002v2))/(vcosφ) φ ′ = ( − p1)/(v2)
$y′ = tan⁡ϕ v′ = -(p1sin⁡ϕ+0.00002⁢v2) vcos⁡ϕ ϕ ′ = -p1v2$
on the range 0 < x < p2$0, with boundary conditions
 y = 0, v = 500, φ = 0.5 at  x = 0, y = 0, v = 450, φ = p3 at  x = p2.
We write y = y(1)$y={\mathbf{y}}\left(1\right)$, v = y(2)$v={\mathbf{y}}\left(2\right)$, φ = y(3)$\varphi ={\mathbf{y}}\left(3\right)$. We estimate p1 = PARAM(1) = 32${p}_{1}=\mathrm{PARAM}\left(1\right)=32$, p2 = PARAM(2) = 6000${p}_{2}=\mathrm{PARAM}\left(2\right)=6000$ and p3 = PARAM(3) = 0.54${p}_{3}=\mathrm{PARAM}\left(3\right)=0.54$ (though this last estimate is not important).
```function nag_ode_bvp_shoot_genpar_example
% Set up initial values for first case.
n = 2;
n1 = 2;
p = [0.2;0];
pe = [1e-05;0.001];
e = [0.0001;0.0001];
m1 = 6;
marray = zeros(1,m1);

fprintf('nag_ode_bvp_shoot_genpar example program results \n');

fprintf('\nCase 1 \n\n');

[pOut1, solnOut1, w1, ifail] = nag_ode_bvp_shoot_genpar(p, pe, e, int64(m1), @fcn1,...
@bc1, @range1, 'n1', int64(n1), 'n', int64(n));
if (ifail ~= 0)
% Problems in integration, or incorrect parameters.  Print
% message and exit.
if (ifail < 5 && ifail ~= 1)
for i = 1:n
fprintf('W(1,2) = %2.2f  ', w(1,2));
fprintf('    W(.,1) = %2.2f',w(i,1));
fprintf('\n');
end
end
error('Warning: nag_ode_bvp_shoot_genpar returned with ifail = %d ',ifail);
else
% Output results.
fprintf('Final parameters \n')
fprintf('    %1.3e   ',pOut1);
fprintf('\n\n Final solution \n')
fprintf('X-value     Components of solution \n')

[x0, x1] = range1(p);
h = (x1-x0)/double(m1-1);
for i = 1:m1; j = 1:n;
m = x0 + double(i-1)*h;
fprintf(' %6.3f        ', m);
fprintf('%2.4f    ',solnOut1(j,i));
fprintf('\n');
% Save results for plotting.
marray(i) = m;
end
end

fprintf('\nCase 2 \n\n');
% Set up initial values for second case.
n = 3;
n1 = 3;
p = [32;6000;0.54];
pe = [1e-05;1e-4;1e-4];
e = [1e-2;1e-2;1e-2];
m1 = 6;
qarray = zeros(1,m1);

[pOut2, solnOut2, w2, ifail] = nag_ode_bvp_shoot_genpar(p, pe, e, int64(m1), @fcn2,...
@bc2, @range2, 'n1', int64(n1), 'n', int64(n));
if (ifail ~= 0)
% Problems in integration, or incorrect parameters.  Print
% message and exit.
if (ifail < 5 && ifail ~= 1)
for i = 1:n
fprintf('W(1,2) = %2.2f  ', w(1,2));
fprintf('    W(.,1) = %2.2f\n',w(i,1));
end
end
error('Warning: nag_ode_bvp_shoot_genpar returned with ifail = %d ',ifail);
else
% Output results.
fprintf('Final parameters \n')
fprintf('    %1.3e   ',pOut2);
fprintf('\n\n Final solution \n')
fprintf(' X-value        Components of solution \n')

[x0, x1] = range2(pOut2);
h = (x1-x0)/double(m1-1);
for i = 1:m1; j = 1:n;
q = x0 + double(i-1)*h;
fprintf('%6.0f     ', q);
fprintf('%10.4f  ',solnOut2(j,i));
fprintf('\n');
% Save results for plotting.
qarray(i) = q;
end
end

% Plot results.
fig = figure('Number', 'off');
display_plot(marray, solnOut1, 'Derivative');

fig = figure('Number', 'off');
display_plot(marray, solnOut2, 'Angle');

function [g1, g2] = bc1(p)
% Evaluate boundary conditions.
z=0.1;
g1(1) = 0.1+p(1)*sqrt(z)*0.1+0.01*z;
g1(2) = p(1)*0.05/sqrt(z) + 0.01;
g2(1) = 1/6;
g2(2) = p(2);
function [g1, g2] = bc2(p)
% Evaluate boundary conditions.
g1(1) = 0.0;
g1(2) = 500.0;
g1(3) = 0.5;
g2(1) = 0.0;
g2(2) = 450.0;
g2(3) = p(3);
function f = fcn1(x, y, p)
% Evaluate derivatives.
f = zeros(2,1);
f(1) = y(2);
f(2) = (y(1)^3-y(2))/2/x;
function f = fcn2(x, y, p)
% Evaluate derivatives.
f = zeros(3,1);
f(1) = tan(y(3));
f(2) = -p(1)*tan(y(3))/y(2) - 0.00002*y(2)/cos(y(3));
f(3) = -p(1)/y(2)^2;
function [x0, x1] = range1(p)
% Evaluate boundary points.
x0 = 0.1;
x1 = 16.0;
function [x0, x1] = range2(p)
% Evaluate boundary points.
x0 = 0.0;
x1 = p(2);
function display_plot(x, y, ylabelString)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Choose the type of plot based on the y axis label.
if strncmp(ylabelString, 'Derivative', 8)
% Plot the results as a linear graph (2 y axes).
[haxes, hline1, hline2] = plotyy(x, y(1,:), x, y(2,:));
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [0.1 0.18]);
set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0.1 0.12 0.14 0.16 0.18]);
set(haxes(2), 'YLim', [0 0.02]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [0 0.005 0.01 0.015 0.02]);
for iaxis = 1:2
% These properties must be the same for both sets of axes.
set(haxes(iaxis), 'XLim', [0 16]);
set(haxes(iaxis), 'XTick', [0 4 8 12 16]);
set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Set the title.
title(['Parameterised Two-point BVP using Initial Value ', ...
'Techniques & Newton Iteration'], titleFmt{:});
% Label the axes.
xlabel('x', labFmt{:});
ylabel(haxes(1),'Solution', labFmt{:});
ylabel(haxes(2),ylabelString, labFmt{:});
legend('solution y(x)','derivative y''(x)','Location','Best');
% Set some features of the two lines.
set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-');
set(hline2, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '--');
else
% Plot the results as a linear graph (2 y axes).
[haxes, hline1, hline2] = plotyy(x, y(1,:), x, y(3,:));
hold on
hline3 = plot(x, y(2,:));
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [-1000 1000]);
set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
set(haxes(1), 'YTick', [-1000 -500 0 500 1000]);
set(haxes(2), 'YLim', [-1 1]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [-1 -0.5 0 0.5 1]);
for iaxis = 1:2
% These properties must be the same for both sets of axes.
set(haxes(iaxis), 'XLim', [0 16]);
set(haxes(iaxis), 'XTick', [0 4 8 12 16]);
set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Set the title.
title(['Find Gravitational Constant & Range from Projectile ', ...
'Terminal Velocity'], titleFmt{:});
% Label the axes.
xlabel('x', labFmt{:});
ylabel(haxes(1),'Height and Velocity', labFmt{:});
ylabel(haxes(2),ylabelString, labFmt{:});
legend('height','velocity','angle','Location','Best');
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-');
set(hline2, 'Linewidth', 0.5, 'Marker', '*', 'LineStyle', ':');
set(hline3, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '--');
end
```
```
nag_ode_bvp_shoot_genpar example program results

Case 1

Final parameters
4.629e-02       3.494e-03

Final solution
X-value     Components of solution
0.100        0.1025    0.0173
3.280        0.1217    0.0042
6.460        0.1338    0.0036
9.640        0.1449    0.0034
12.820        0.1557    0.0034
16.000        0.1667    0.0035

Case 2

Final parameters
3.239e+01       5.962e+03       -5.353e-01

Final solution
X-value        Components of solution
0         0.0000    500.0000      0.5000
1192       529.6471    451.5554      0.3280
2385       807.2140    420.3050      0.1231
3577       820.3972    409.4254     -0.1032
4769       556.1322    419.9890     -0.3296
5962        -0.0003    450.0000     -0.5353

```
```function d02hb_example
% Set up initial values for first case.
n = 2;
n1 = 2;
p = [0.2;0];
pe = [1e-05;0.001];
e = [0.0001;0.0001];
m1 = 6;
marray = zeros(1,m1);

fprintf('d02hb example program results \n');

fprintf('\nCase 1 \n\n');

[pOut1, solnOut1, w1, ifail] = d02hb(p, pe, e, int64(m1), @fcn1,...
@bc1, @range1, 'n1', int64(n1), 'n', int64(n));
if (ifail ~= 0)
% Problems in integration, or incorrect parameters.  Print
% message and exit.
if (ifail < 5 && ifail ~= 1)
for i = 1:n
fprintf('W(1,2) = %2.2f  ', w(1,2));
fprintf('    W(.,1) = %2.2f',w(i,1));
fprintf('\n');
end
end
error('Warning: d02hb returned with ifail = %d ',ifail);
else
% Output results.
fprintf('Final parameters \n')
fprintf('    %1.3e   ',pOut1);
fprintf('\n\n Final solution \n')
fprintf('X-value     Components of solution \n')

[x0, x1] = range1(p);
h = (x1-x0)/double(m1-1);
for i = 1:m1; j = 1:n;
m = x0 + double(i-1)*h;
fprintf(' %6.3f        ', m);
fprintf('%2.4f    ',solnOut1(j,i));
fprintf('\n');
% Save results for plotting.
marray(i) = m;
end
end

fprintf('\nCase 2 \n\n');
% Set up initial values for second case.
n = 3;
n1 = 3;
p = [32;6000;0.54];
pe = [1e-05;1e-4;1e-4];
e = [1e-2;1e-2;1e-2];
m1 = 6;
qarray = zeros(1,m1);

[pOut2, solnOut2, w2, ifail] = d02hb(p, pe, e, int64(m1), @fcn2,...
@bc2, @range2, 'n1', int64(n1), 'n', int64(n));
if (ifail ~= 0)
% Problems in integration, or incorrect parameters.  Print
% message and exit.
if (ifail < 5 && ifail ~= 1)
for i = 1:n
fprintf('W(1,2) = %2.2f  ', w(1,2));
fprintf('    W(.,1) = %2.2f\n',w(i,1));
end
end
error('Warning: d02hb returned with ifail = %d ',ifail);
else
% Output results.
fprintf('Final parameters \n')
fprintf('    %1.3e   ',pOut2);
fprintf('\n\n Final solution \n')
fprintf(' X-value        Components of solution \n')

[x0, x1] = range2(pOut2);
h = (x1-x0)/double(m1-1);
for i = 1:m1; j = 1:n;
q = x0 + double(i-1)*h;
fprintf('%6.0f     ', q);
fprintf('%10.4f  ',solnOut2(j,i));
fprintf('\n');
% Save results for plotting.
qarray(i) = q;
end
end

% Plot results.
fig = figure('Number', 'off');
display_plot(marray, solnOut1, 'Derivative');

fig = figure('Number', 'off');
display_plot(marray, solnOut2, 'Angle');

function [g1, g2] = bc1(p)
% Evaluate boundary conditions.
z=0.1;
g1(1) = 0.1+p(1)*sqrt(z)*0.1+0.01*z;
g1(2) = p(1)*0.05/sqrt(z) + 0.01;
g2(1) = 1/6;
g2(2) = p(2);
function [g1, g2] = bc2(p)
% Evaluate boundary conditions.
g1(1) = 0.0;
g1(2) = 500.0;
g1(3) = 0.5;
g2(1) = 0.0;
g2(2) = 450.0;
g2(3) = p(3);
function f = fcn1(x, y, p)
% Evaluate derivatives.
f = zeros(2,1);
f(1) = y(2);
f(2) = (y(1)^3-y(2))/2/x;
function f = fcn2(x, y, p)
% Evaluate derivatives.
f = zeros(3,1);
f(1) = tan(y(3));
f(2) = -p(1)*tan(y(3))/y(2) - 0.00002*y(2)/cos(y(3));
f(3) = -p(1)/y(2)^2;
function [x0, x1] = range1(p)
% Evaluate boundary points.
x0 = 0.1;
x1 = 16.0;
function [x0, x1] = range2(p)
% Evaluate boundary points.
x0 = 0.0;
x1 = p(2);
function display_plot(x, y, ylabelString)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Choose the type of plot based on the y axis label.
if strncmp(ylabelString, 'Derivative', 8)
% Plot the results as a linear graph (2 y axes).
[haxes, hline1, hline2] = plotyy(x, y(1,:), x, y(2,:));
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [0.1 0.18]);
set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0.1 0.12 0.14 0.16 0.18]);
set(haxes(2), 'YLim', [0 0.02]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [0 0.005 0.01 0.015 0.02]);
for iaxis = 1:2
% These properties must be the same for both sets of axes.
set(haxes(iaxis), 'XLim', [0 16]);
set(haxes(iaxis), 'XTick', [0 4 8 12 16]);
set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Set the title.
title(['Parameterised Two-point BVP using Initial Value ', ...
'Techniques & Newton Iteration'], titleFmt{:});
% Label the axes.
xlabel('x', labFmt{:});
ylabel(haxes(1),'Solution', labFmt{:});
ylabel(haxes(2),ylabelString, labFmt{:});
legend('solution y(x)','derivative y''(x)','Location','Best');
% Set some features of the two lines.
set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-');
set(hline2, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '--');
else
% Plot the results as a linear graph (2 y axes).
[haxes, hline1, hline2] = plotyy(x, y(1,:), x, y(3,:));
hold on
hline3 = plot(x, y(2,:));
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [-1000 1000]);
set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
set(haxes(1), 'YTick', [-1000 -500 0 500 1000]);
set(haxes(2), 'YLim', [-1 1]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [-1 -0.5 0 0.5 1]);
for iaxis = 1:2
% These properties must be the same for both sets of axes.
set(haxes(iaxis), 'XLim', [0 16]);
set(haxes(iaxis), 'XTick', [0 4 8 12 16]);
set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Set the title.
title(['Find Gravitational Constant & Range from Projectile ', ...
'Terminal Velocity'], titleFmt{:});
% Label the axes.
xlabel('x', labFmt{:});
ylabel(haxes(1),'Height and Velocity', labFmt{:});
ylabel(haxes(2),ylabelString, labFmt{:});
legend('height','velocity','angle','Location','Best');
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-');
set(hline2, 'Linewidth', 0.5, 'Marker', '*', 'LineStyle', ':');
set(hline3, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '--');
end
```
```
d02hb example program results

Case 1

Final parameters
4.629e-02       3.494e-03

Final solution
X-value     Components of solution
0.100        0.1025    0.0173
3.280        0.1217    0.0042
6.460        0.1338    0.0036
9.640        0.1449    0.0034
12.820        0.1557    0.0034
16.000        0.1667    0.0035

Case 2

Final parameters
3.239e+01       5.962e+03       -5.353e-01

Final solution
X-value        Components of solution
0         0.0000    500.0000      0.5000
1192       529.6471    451.5554      0.3280
2385       807.2140    420.3050      0.1231
3577       820.3972    409.4254     -0.1032
4769       556.1322    419.9890     -0.3296
5962        -0.0003    450.0000     -0.5353

```