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_sum_invlaplace_weeks (c06lb)

## Purpose

nag_sum_invlaplace_weeks (c06lb) computes the inverse Laplace transform f(t) $f\left(t\right)$ of a user-supplied function F(s) $F\left(s\right)$, defined for complex s$s$. The function uses a modification of Weeks' method which is suitable when f(t) $f\left(t\right)$ has continuous derivatives of all orders. The function returns the coefficients of an expansion which approximates f(t) $f\left(t\right)$ and can be evaluated for given values of t$t$ by subsequent calls of nag_sum_invlaplace_weeks_eval (c06lc).

## Syntax

[sigma, b, m, acoef, errvec, ifail] = c06lb(f, sigma0, sigma, b, epstol, 'mmax', mmax)
[sigma, b, m, acoef, errvec, ifail] = nag_sum_invlaplace_weeks(f, sigma0, sigma, b, epstol, 'mmax', mmax)

## Description

Given a function f (t) $f\left(t\right)$ of a real variable t$t$, its Laplace transform F(s) $F\left(s\right)$ is a function of a complex variable s$s$, defined by
 ∞ F(s) = ∫ e − stf(t)dt,  Re(s) > σ0. 0
$F (s) = ∫0∞ e-st f (t) dt , Re(s) > σ0 .$
Then f(t) $f\left(t\right)$ is the inverse Laplace transform of F(s) $F\left(s\right)$. The value σ0 ${\sigma }_{0}$ is referred to as the abscissa of convergence of the Laplace transform; it is the rightmost real part of the singularities of F(s) $F\left(s\right)$.
nag_sum_invlaplace_weeks (c06lb), along with its companion nag_sum_invlaplace_weeks_eval (c06lc), attempts to solve the following problem:
• given a function F(s) $F\left(s\right)$, compute values of its inverse Laplace transform f(t) $f\left(t\right)$ for specified values of t$t$.
The method is a modification of Weeks' method (see Garbow et al. (1988a)), which approximates f(t) $f\left(t\right)$ by a truncated Laguerre expansion:
 m − 1 f̃(t) = eσt ∑ aie − bt / 2 Li(bt),  σ > σ0,  b > 0 i = 0
$f~ (t) = eσt ∑ i=0 m-1 ai e -bt / 2 Li (bt) , σ > σ0 , b > 0$
where Li (x) ${L}_{i}\left(x\right)$ is the Laguerre polynomial of degree i$i$. This function computes the coefficients ai ${a}_{i}$ of the above Laguerre expansion; the expansion can then be evaluated for specified t$t$ by calling nag_sum_invlaplace_weeks_eval (c06lc). You must supply the value of σ0 ${\sigma }_{0}$, and also suitable values for σ$\sigma$ and b$b$: see Section [Further Comments] for guidance.
The method is only suitable when f(t) $f\left(t\right)$ has continuous derivatives of all orders. For such functions the approximation (t) $\stackrel{~}{f}\left(t\right)$ is usually good and inexpensive. The function will fail with an error exit if the method is not suitable for the supplied function F(s) $F\left(s\right)$.
The function is designed to satisfy an accuracy criterion of the form:
 |( f(t) − f̃(t) )/(eσt)| < εtol,   for all ​t $| f(t)- f~(t) e σt | < ε tol , for all ​t$
where εtol ${\epsilon }_{\mathit{tol}}$ is a user-supplied bound. The error measure on the left-hand side is referred to as the pseudo-relative error, or pseudo-error for short. Note that if σ > 0 $\sigma >0$ and t$t$ is large, the absolute error in (t) $\stackrel{~}{f}\left(t\right)$ may be very large.
nag_sum_invlaplace_weeks (c06lb) is derived from the function MODUL1 in Garbow et al. (1988a).

## References

Garbow B S, Giunta G, Lyness J N and Murli A (1988a) Software for an implementation of Weeks' method for the inverse laplace transform problem ACM Trans. Math. Software 14 163–170
Garbow B S, Giunta G, Lyness J N and Murli A (1988b) Algorithm 662: A Fortran software package for the numerical inversion of the Laplace transform based on Weeks' method ACM Trans. Math. Software 14 171–176

## Parameters

### Compulsory Input Parameters

1:     f – function handle or string containing name of m-file
f must return the value of the Laplace transform function F(s)$F\left(s\right)$ for a given complex value of s$s$.
[result] = f(s)

Input Parameters

1:     s – complex scalar
The value of s$s$ for which F(s)$F\left(s\right)$ must be evaluated. The real part of s is greater than σ0${\sigma }_{0}$.

Output Parameters

1:     result – complex scalar
The result of the function.
2:     sigma0 – double scalar
The abscissa of convergence of the Laplace transform, σ0${\sigma }_{0}$.
3:     sigma – double scalar
The parameter σ$\sigma$ of the Laguerre expansion. If on entry sigmaσ0${\mathbf{sigma}}\le {\sigma }_{0}$, sigma is reset to σ0 + 0.7${\sigma }_{0}+0.7$.
4:     b – double scalar
The parameter b$b$ of the Laguerre expansion. If on entry b < 2 (σσ0)${\mathbf{b}}<2\left(\sigma -{\sigma }_{0}\right)$, b is reset to 2.5(σσ0)$2.5\left(\sigma -{\sigma }_{0}\right)$.
5:     epstol – double scalar
The required relative pseudo-accuracy, that is, an upper bound on |f(t)(t)| eσt$|f\left(t\right)-\stackrel{~}{f}\left(t\right)|{e}^{-\sigma t}$.

### Optional Input Parameters

1:     mmax – int64int32nag_int scalar
An upper bound on the number of Laguerre expansion coefficients to be computed. The number of coefficients actually computed is always a power of 2$2$, so mmax should be a power of 2$2$; if mmax is not a power of 2$2$ then the maximum number of coefficients calculated will be the largest power of 2$2$ less than mmax.
Default: 1024$1024$
Constraint: mmax8${\mathbf{mmax}}\ge 8$.

None.

### Output Parameters

1:     sigma – double scalar
The value actually used for σ$\sigma$, as just described.
2:     b – double scalar
The value actually used for b$b$, as just described.
3:     m – int64int32nag_int scalar
The number of Laguerre expansion coefficients actually computed. The number of calls to f is m / 2 + 2${\mathbf{m}}/2+2$.
4:     acoef(mmax) – double array
The first m elements contain the computed Laguerre expansion coefficients, ai${a}_{i}$.
5:     errvec(8$8$) – double array
An 8$8$-component vector of diagnostic information.
errvec(1)${\mathbf{errvec}}\left(1\right)$
Overall estimate of the pseudo-error
|f(t)(t)| eσt = errvec(2) + errvec(3) + errvec(4)$|f\left(t\right)-\stackrel{~}{f}\left(t\right)|{e}^{-\sigma t}={\mathbf{errvec}}\left(2\right)+{\mathbf{errvec}}\left(3\right)+{\mathbf{errvec}}\left(4\right)$.
errvec(2)${\mathbf{errvec}}\left(2\right)$
Estimate of the discretization pseudo-error.
errvec(3)${\mathbf{errvec}}\left(3\right)$
Estimate of the truncation pseudo-error.
errvec(4)${\mathbf{errvec}}\left(4\right)$
Estimate of the condition pseudo-error on the basis of minimal noise levels in function values.
errvec(5)${\mathbf{errvec}}\left(5\right)$
K$K$, coefficient of a heuristic decay function for the expansion coefficients.
errvec(6)${\mathbf{errvec}}\left(6\right)$
R$R$, base of the decay function for the expansion coefficients.
errvec(7)${\mathbf{errvec}}\left(7\right)$
Logarithm of the largest expansion coefficient.
errvec(8)${\mathbf{errvec}}\left(8\right)$
Logarithm of the smallest nonzero expansion coefficient.
The values K$K$ and R$R$ returned in errvec(5)${\mathbf{errvec}}\left(5\right)$ and errvec(6)${\mathbf{errvec}}\left(6\right)$ define a decay function KRi$K{R}^{-i}$ constructed by the function for the purposes of error estimation. It satisfies
 |ai| < K R − i ,   ​ i = 1, 2, … , m . $|ai| < K R -i , ​ i= 1, 2, …, m .$
6:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Note: nag_sum_invlaplace_weeks (c06lb) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

ifail = 1${\mathbf{ifail}}=1$
 On entry, mmax < 8${\mathbf{mmax}}<8$.
W ifail = 2${\mathbf{ifail}}=2$
The estimated pseudo-error bounds are slightly larger than epstol. Note, however, that the actual errors in the final results may be smaller than epstol as bounds independent of the value of t$t$ are pessimistic.
W ifail = 3${\mathbf{ifail}}=3$
Computation was terminated early because the estimate of rounding error was greater than epstol. Increasing epstol may help.
ifail = 4${\mathbf{ifail}}=4$
The decay rate of the coefficients is too small. Increasing mmax may help.
ifail = 5${\mathbf{ifail}}=5$
The decay rate of the coefficients is too small. In addition the rounding error is such that the required accuracy cannot be obtained. Increasing mmax or epstol may help.
W ifail = 6${\mathbf{ifail}}=6$
The behaviour of the coefficients does not enable reasonable prediction of error bounds. Check the value of sigma0. In this case, errvec(i) ${\mathbf{errvec}}\left(\mathit{i}\right)$ is set to 1.0 $-1.0$, for i = 1,2,,5$\mathit{i}=1,2,\dots ,5$.
When ${\mathbf{ifail}}\ge {\mathbf{3}}$, changing sigma or b may help. If not, the method should be abandoned.

## Accuracy

The error estimate returned in errvec(1) ${\mathbf{errvec}}\left(1\right)$ has been found in practice to be a highly reliable bound on the pseudo-error |f(t)(t)| eσt $|f\left(t\right)-\stackrel{~}{f}\left(t\right)|{e}^{-\sigma t}$.

### The Role of σ0

Nearly all techniques for inversion of the Laplace transform require you to supply the value of σ0 ${\sigma }_{0}$, the convergence abscissa, or else an upper bound on σ0 ${\sigma }_{0}$. For this function, one of the reasons for having to supply σ0 ${\sigma }_{0}$ is that the parameter σ$\sigma$ must be greater than σ0 ${\sigma }_{0}$; otherwise the series for (t) $\stackrel{~}{f}\left(t\right)$ will not converge.
If you do not know the value of σ0 ${\sigma }_{0}$, you must be prepared for significant preliminary effort, either in experimenting with the method and obtaining chaotic results, or in attempting to locate the rightmost singularity of F(s) $F\left(s\right)$.
The value of σ0 ${\sigma }_{0}$ is also relevant in defining a natural accuracy criterion. For large t$t$, f(t) $f\left(t\right)$ is of uniform numerical order k eσ0t $k{e}^{{\sigma }_{0}t}$, so a natural measure of relative accuracy of the approximation (t) $\stackrel{~}{f}\left(t\right)$ is:
 εnat (t) = ( f̃ (t) − f (t) ) / eσ0t . $εnat (t) = ( f~ (t) - f (t) ) / e σ0t .$
nag_sum_invlaplace_weeks (c06lb) uses the supplied value of σ0 ${\sigma }_{0}$ only in determining the values of σ$\sigma$ and b$b$ (see Sections [Choice of ] and [Choice of ]); thereafter it bases its computation entirely on σ$\sigma$ and b$b$.

### Choice of σ

Even when the value of σ0 ${\sigma }_{0}$ is known, choosing a value for σ$\sigma$ is not easy. Briefly, the series for (t) $\stackrel{~}{f}\left(t\right)$ converges slowly when σσ0 $\sigma -{\sigma }_{0}$ is small, and faster when σσ0 $\sigma -{\sigma }_{0}$ is larger. However the natural accuracy measure satisfies
 |εnat(t)| < εtol e (σ − σ0) t $| εnat (t) | < εtol e ( σ - σ0 ) t$
and this degrades exponentially with t$t$, the exponential constant being σσ0 $\sigma -{\sigma }_{0}$.
Hence, if you require meaningful results over a large range of values of t$t$, you should choose σσ0 $\sigma -{\sigma }_{0}$ small, in which case the series for (t) $\stackrel{~}{f}\left(t\right)$ converges slowly; while for a smaller range of values of t$t$, you can allow σσ0 $\sigma -{\sigma }_{0}$ to be larger and obtain faster convergence.
The default value for σ$\sigma$ used by nag_sum_invlaplace_weeks (c06lb) is σ0 + 0.7 ${\sigma }_{0}+0.7$. There is no theoretical justification for this.

### Choice of b

The simplest advice for choosing b$b$ is to set b / 2 σ σ0 $b/2\ge \sigma -{\sigma }_{0}$. The default value used by the function is 2.5 (σσ0) $2.5\left(\sigma -{\sigma }_{0}\right)$. A more refined choice is to set
 b / 2 ≥ min |σ − sj| j
$b/2 ≥ minj |σ-sj|$
where sj ${s}_{j}$ are the singularities of F(s) $F\left(s\right)$.

## Example

This example computes values of the inverse Laplace transform of the function
 F(s) = 3/(s2 − 9) . $F(s) = 3 s2-9 .$
 f(t) = sinh3t . $f(t) = sinh⁡3t .$
The program first calls nag_sum_invlaplace_weeks (c06lb) to compute the coefficients of the Laguerre expansion, and then calls nag_sum_invlaplace_weeks_eval (c06lc) to evaluate the expansion at t = 0$t=0$, 1$1$, 2$2$, 3$3$, 4$4$, 5$5$.
```function nag_sum_invlaplace_weeks_example
% Initialize variables and arrays.
sigma0 = 3;
epstol = 1e-5;
b = 0;
sigma = 0;
n = 5;
earray = zeros(1, n+1);
jarray = zeros(1, n+1);
farray = zeros(1, n+1);
parray = zeros(1, n+1);

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

[sigmaOut, bOut, m, acoef, errvec, ifail] = ...
nag_sum_invlaplace_weeks(@f, sigma0, sigma, b, epstol);
if ifail ~= 0
% Convergence problems.  Print message and exit.
error('Warning: nag_sum_invlaplace_weeks returned with ifail = %1d ',ifail);
end

% Prepare to output results.
disp(['No. of coefficients returned by nag_sum_invlaplace_weeks = ',num2str(m)]);
disp('  ');
disp('      Computed       Exact         Pseudo');
disp('t       f(t)          f(t)          error');
% Evaluate inverse transform for different values of t.  We use nag_sum_invlaplace_weeks_eval,
% which calculates the transform from the coefficients given by nag_sum_invlaplace_weeks.
for j = 0:5
t = double(j);

[finv, ifail] = nag_sum_invlaplace_weeks_eval(t, sigmaOut, bOut, acoef, errvec, 'm', m);
if ifail ~= 0
% Approximation is too large or too small.  Print message and exit.
error('Warning: nag_sum_invlaplace_weeks_eval returned with ifail = %1d ',ifail);
end
exact = sinh(3.0*t);
pserr = abs(finv-exact)/exp(sigmaOut*t);
fprintf('%d   %10.4d   %11.4d     %8.4d\n', t, finv, exact, pserr);
% Create arrays to be used for plotting.
jarray(j+1) = t;
farray(j+1) = finv;
parray(j+1) = pserr;
end

% Plot results.
fig = figure('Number', 'off');
display_plot(jarray, farray, parray)

function [f] = f(s)
% Evaluate the Laplace transform function.
f=3.0/(s^2-9.0);
if isreal(f)
f=complex(f);
end
function plot(jarray, farray, parray)
% 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.
% Use a log plot for both curves.
[haxes, hline1, hline2] = plotyy(jarray, farray, jarray, parray,...
'semilogy','semilogy');
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [1.0e-10 1.0e10]);
set(haxes(1), 'YMinorTick', 'on');
set(haxes(1), 'YTick', [1.0e-10 1.0e-5 1.0 1.0e5 1.0e10]);
set(haxes(2), 'YLim', [1.0e-10 1.0e-8]);
%set(haxes(2), 'FontSize', 11);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [1e-10 1e-9 1e-8]);
for iaxis = 1:2
% These properties must be the same for both sets of axes.
set(haxes(iaxis), 'XLim', [0 5]);
set(haxes(iaxis), 'XTick', [0 1 2 3 4 5]);
set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Set the title.
title('Inverse Laplace Transform of 3/(s^2-9)', titleFmt{:});
% Label the x axis.
xlabel('t', labFmt{:});
% Label the left & right y axes.
ylabel(haxes(1),'f(t)', labFmt{:});
ylabel(haxes(2),'Pseudo Error', labFmt{:});
% Label the curves.
legend('f(t)','Pseudo Error','location','North')
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-');
set(hline2, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '-');
function display_plot(jarray, farray, parray)
% 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.
% Use a log plot for both curves.
[haxes, hline1, hline2] = plotyy(jarray, farray, jarray, parray,...
'semilogy','semilogy');
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [1.0e-10 1.0e10]);
set(haxes(1), 'YMinorTick', 'on');
set(haxes(1), 'YTick', [1.0e-10 1.0e-5 1.0 1.0e5 1.0e10]);
set(haxes(2), 'YLim', [1.0e-10 1.0e-8]);
%set(haxes(2), 'FontSize', 11);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [1e-10 1e-9 1e-8]);
for iaxis = 1:2
% These properties must be the same for both sets of axes.
set(haxes(iaxis), 'XLim', [0 5]);
set(haxes(iaxis), 'XTick', [0 1 2 3 4 5]);
set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Set the title.
title('Inverse Laplace Transform of 3/(s^2-9)', titleFmt{:});
% Label the x axis.
xlabel('t', labFmt{:});
% Label the left & right y axes.
ylabel(haxes(1),'f(t)', labFmt{:});
ylabel(haxes(2),'Pseudo Error', labFmt{:});
% Label the curves.
legend('f(t)','Pseudo Error','location','North')
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-');
set(hline2, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '-');
```
```
nag_sum_invlaplace_weeks example program results

No. of coefficients returned by nag_sum_invlaplace_weeks = 64

Computed       Exact         Pseudo
t       f(t)          f(t)          error
0   1.5129e-09          0000     1.5129e-09
1   1.0018e+01    1.0018e+01     1.7394e-09
2   2.0171e+02    2.0171e+02     1.2471e-10
3   4.0515e+03    4.0515e+03     9.7722e-10
4   8.1377e+04    8.1377e+04     3.0221e-10
5   1.6345e+06    1.6345e+06     1.6991e-09

```
```function c06lb_example
% Initialize variables and arrays.
sigma0 = 3;
epstol = 1e-5;
b = 0;
sigma = 0;
n = 5;
earray = zeros(1, n+1);
jarray = zeros(1, n+1);
farray = zeros(1, n+1);
parray = zeros(1, n+1);

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

[sigmaOut, bOut, m, acoef, errvec, ifail] = ...
c06lb(@f, sigma0, sigma, b, epstol);
if ifail ~= 0
% Convergence problems.  Print message and exit.
error('Warning: c06lb returned with ifail = %1d ',ifail);
end

% Prepare to output results.
disp(['No. of coefficients returned by c06lb = ',num2str(m)]);
disp('  ');
disp('      Computed       Exact         Pseudo');
disp('t       f(t)          f(t)          error');
% Evaluate inverse transform for different values of t.  We use c06lc,
% which calculates the transform from the coefficients given by c06lb.
for j = 0:5
t = double(j);

[finv, ifail] = c06lc(t, sigmaOut, bOut, acoef, errvec, 'm', m);
if ifail ~= 0
% Approximation is too large or too small.  Print message and exit.
error('Warning: c06lc returned with ifail = %1d ',ifail);
end
exact = sinh(3.0*t);
pserr = abs(finv-exact)/exp(sigmaOut*t);
fprintf('%d   %10.4d   %11.4d     %8.4d\n', t, finv, exact, pserr);
% Create arrays to be used for plotting.
jarray(j+1) = t;
farray(j+1) = finv;
parray(j+1) = pserr;
end

% Plot results.
fig = figure('Number', 'off');
display_plot(jarray, farray, parray)

function [f] = f(s)
% Evaluate the Laplace transform function.
f=3.0/(s^2-9.0);
if isreal(f)
f=complex(f);
end
function plot(jarray, farray, parray)
% 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.
% Use a log plot for both curves.
[haxes, hline1, hline2] = plotyy(jarray, farray, jarray, parray,...
'semilogy','semilogy');
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [1.0e-10 1.0e10]);
set(haxes(1), 'YMinorTick', 'on');
set(haxes(1), 'YTick', [1.0e-10 1.0e-5 1.0 1.0e5 1.0e10]);
set(haxes(2), 'YLim', [1.0e-10 1.0e-8]);
%set(haxes(2), 'FontSize', 11);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [1e-10 1e-9 1e-8]);
for iaxis = 1:2
% These properties must be the same for both sets of axes.
set(haxes(iaxis), 'XLim', [0 5]);
set(haxes(iaxis), 'XTick', [0 1 2 3 4 5]);
set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Set the title.
title('Inverse Laplace Transform of 3/(s^2-9)', titleFmt{:});
% Label the x axis.
xlabel('t', labFmt{:});
% Label the left & right y axes.
ylabel(haxes(1),'f(t)', labFmt{:});
ylabel(haxes(2),'Pseudo Error', labFmt{:});
% Label the curves.
legend('f(t)','Pseudo Error','location','North')
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-');
set(hline2, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '-');
function display_plot(jarray, farray, parray)
% 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.
% Use a log plot for both curves.
[haxes, hline1, hline2] = plotyy(jarray, farray, jarray, parray,...
'semilogy','semilogy');
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [1.0e-10 1.0e10]);
set(haxes(1), 'YMinorTick', 'on');
set(haxes(1), 'YTick', [1.0e-10 1.0e-5 1.0 1.0e5 1.0e10]);
set(haxes(2), 'YLim', [1.0e-10 1.0e-8]);
%set(haxes(2), 'FontSize', 11);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [1e-10 1e-9 1e-8]);
for iaxis = 1:2
% These properties must be the same for both sets of axes.
set(haxes(iaxis), 'XLim', [0 5]);
set(haxes(iaxis), 'XTick', [0 1 2 3 4 5]);
set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Set the title.
title('Inverse Laplace Transform of 3/(s^2-9)', titleFmt{:});
% Label the x axis.
xlabel('t', labFmt{:});
% Label the left & right y axes.
ylabel(haxes(1),'f(t)', labFmt{:});
ylabel(haxes(2),'Pseudo Error', labFmt{:});
% Label the curves.
legend('f(t)','Pseudo Error','location','North')
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-');
set(hline2, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '-');
```
```
c06lb example program results

No. of coefficients returned by c06lb = 64

Computed       Exact         Pseudo
t       f(t)          f(t)          error
0   1.5129e-09          0000     1.5129e-09
1   1.0018e+01    1.0018e+01     1.7394e-09
2   2.0171e+02    2.0171e+02     1.2471e-10
3   4.0515e+03    4.0515e+03     9.7722e-10
4   8.1377e+04    8.1377e+04     3.0221e-10
5   1.6345e+06    1.6345e+06     1.6991e-09

```