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_dae_dassl_gen (d02ne)

## Purpose

nag_ode_dae_dassl_gen (d02ne) is a function for integrating stiff systems of implicit ordinary differential equations coupled with algebraic equations.

## Syntax

[t, y, ydot, rtol, atol, itask, icom, com, user, ifail] = d02ne(t, tout, y, ydot, rtol, atol, itask, res, jac, icom, com, 'neq', neq, 'user', user)
[t, y, ydot, rtol, atol, itask, icom, com, user, ifail] = nag_ode_dae_dassl_gen(t, tout, y, ydot, rtol, atol, itask, res, jac, icom, com, 'neq', neq, 'user', user)

## Description

nag_ode_dae_dassl_gen (d02ne) is a general purpose function for integrating the initial value problem for a stiff system of implicit ordinary differential equations with coupled algebraic equations written in the form
 F (t,y,y′) = 0 . $F (t,y,y′) = 0 .$
nag_ode_dae_dassl_gen (d02ne) uses the DASSL implementation of the Backward Differentiation Formulae (BDF) of orders one to five to solve a system of the above form for y$y$ (y) and y${y}^{\prime }$ (ydot). Values for y and ydot at the initial time must be given as input. These values must be consistent, (i.e., if t, y, ydot are the given initial values, they must satisfy F(t,y,ydot) = 0$F\left({\mathbf{t}},{\mathbf{y}},{\mathbf{ydot}}\right)=0$). The function solves the system from t = t$t={\mathbf{t}}$ to t = tout$t={\mathbf{tout}}$.
An outline of a typical calling program for nag_ode_dae_dassl_gen (d02ne) is given below. It calls the DASSL implementation of the BDF integrator setup function nag_ode_dae_dassl_setup (d02mw) and the banded matrix setup function nag_ode_dae_dassl_linalg (d02np) (if required), and, if the integration needs to proceed, calls nag_ode_dae_dassl_cont (d02mc) before continuing the integration.
```.
.
.
% Initialize the integrator
[...] = d02mw(...);
% Is the Jacobian matrix banded?
if (banded)
[...] = d02np(...);
end
% Set DT to the required temporal resolution
% Set TEND to the final time
% Call the integrator for each temporal value
while 1
[tout,...] = d02ne(...);
break
end
tout = min(tout + dt, tend);
end
% Print the solution
[...] = d02mc(...);
end
.
.
.
```

None.

## Parameters

### Compulsory Input Parameters

1:     t – double scalar
On initial entry: the initial value of the independent variable, t$t$.
2:     tout – double scalar
The next value of t$t$ at which a computed solution is desired.
On initial entry: tout is used to determine the direction of integration. Integration is permitted in either direction (see also itask).
Constraint: ${\mathbf{tout}}\ne {\mathbf{t}}$.
3:     y(neq) – double array
neq, the dimension of the array, must satisfy the constraint neq1${\mathbf{neq}}\ge 1$.
On initial entry: the vector of initial values of the dependent variables y$y$.
4:     ydot(neq) – double array
neq, the dimension of the array, must satisfy the constraint neq1${\mathbf{neq}}\ge 1$.
On initial entry: ydot must contain approximations to the time derivatives y ' $y\text{'}$ of the vector y$y$ evaluated at the initial value of the independent variable.
5:     rtol( : $:$) – double array
Note: the dimension of the array rtol depends on the value of itol as set in nag_ode_dae_dassl_setup (d02mw); it  must be at least neq${\mathbf{neq}}$ if itol = true${\mathbf{itol}}=\mathbf{true}$ and at least 1$1$ if itol = false${\mathbf{itol}}=\mathbf{false}$.
The relative local error tolerance.
Constraint: rtol(i)0.0${\mathbf{rtol}}\left(\mathit{i}\right)\ge 0.0$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$
where n = neq$n={\mathbf{neq}}$ when itol = true${\mathbf{itol}}=\mathbf{true}$ and n = 1$n=1$ otherwise.
6:     atol( : $:$) – double array
Note: the dimension of the array atol depends on the value of itol as set in nag_ode_dae_dassl_setup (d02mw); it  must be at least neq${\mathbf{neq}}$ if itol = true${\mathbf{itol}}=\mathbf{true}$ and at least 1$1$ if itol = false${\mathbf{itol}}=\mathbf{false}$.
The absolute local error tolerance.
Constraint: atol(i)0.0${\mathbf{atol}}\left(\mathit{i}\right)\ge 0.0$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$
where n = neq$n={\mathbf{neq}}$ when itol = true${\mathbf{itol}}=\mathbf{true}$ and n = 1$n=1$ otherwise.
On initial entry: need not be set.
8:     res – function handle or string containing name of m-file
res must evaluate the residual
 R = F (t,y,y′) . $R = F (t,y,y′) .$
[r, ires, user] = res(neq, t, y, ydot, ires, user)

Input Parameters

1:     neq – int64int32nag_int scalar
The number of differential-algebraic equations being solved.
2:     t – double scalar
t$t$, the current value of the independent variable.
3:     y(neq) – double array
yi${y}_{\mathit{i}}$, for i = 1,2,,neq$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, the current solution component.
4:     ydot(neq) – double array
The derivative of the solution at the current point t$t$.
5:     ires – int64int32nag_int scalar
Is always equal to zero.
6:     user – Any MATLAB object
res is called from nag_ode_dae_dassl_gen (d02ne) with the object supplied to nag_ode_dae_dassl_gen (d02ne).

Output Parameters

1:     r(neq) – double array
r(i)${\mathbf{r}}\left(\mathit{i}\right)$ must contain the i$\mathit{i}$th component of R$R$, for i = 1,2,,neq$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$ where
 R = F (t,y,ydot) . $R = F (t,y,ydot) .$
2:     ires – int64int32nag_int scalar
ires should normally be left unchanged. However, if an illegal value of y is encountered, ires should be set to 1$-1$; nag_ode_dae_dassl_gen (d02ne) will then attempt to resolve the problem so that illegal values of y are not encountered. ires should be set to 2$-2$ if you wish to return control to the calling (sub)routine; this will cause nag_ode_dae_dassl_gen (d02ne) to exit with ${\mathbf{ifail}}={\mathbf{23}}$.
3:     user – Any MATLAB object
9:     jac – function handle or string containing name of m-file
Evaluates the matrix of partial derivatives, J$J$, where
 Jij = (∂Fi)/(∂yj) + cj × (∂Fi)/(∂y′j) ,  i,j = 1,2, … ,neq . $J ij = ∂Fi ∂yj + cj × ∂Fi ∂y′j , i,j=1,2,…,neq .$
If this option is not required, the actual argument for jac must be the string 'd02nez'. (nag_ode_dae_dassl_gen_dummy_jac (d02nez) is included in the NAG Toolbox.) You must indicate to the integrator whether this option is to be used by setting the parameter jceval appropriately in a call to the setup function nag_ode_dae_dassl_setup (d02mw).
[pd, user] = jac(neq, t, y, ydot, pd, cj, user)

Input Parameters

1:     neq – int64int32nag_int scalar
The number of differential-algebraic equations being solved.
2:     t – double scalar
t$t$, the current value of the independent variable.
3:     y(neq) – double array
yi${y}_{\mathit{i}}$, for i = 1,2,,neq$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, the current solution component.
4:     ydot(neq) – double array
The derivative of the solution at the current point t$t$.
5:     pd( : $:$) – double array
pd is preset to zero before the call to jac.
6:     cj – double scalar
cj is a scalar constant which will be defined in nag_ode_dae_dassl_gen (d02ne).
7:     user – Any MATLAB object
jac is called from nag_ode_dae_dassl_gen (d02ne) with the object supplied to nag_ode_dae_dassl_gen (d02ne).

Output Parameters

1:     pd( : $:$) – double array
If the Jacobian is full then pd((ij) × neq + i) = Jij${\mathbf{pd}}\left(\left(\mathit{i}-\mathit{j}\right)×{\mathbf{neq}}+\mathit{i}\right)={J}_{\mathit{i}\mathit{j}}$, for i = 1,2,,neq$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$ and j = 1,2,,neq$\mathit{j}=1,2,\dots ,{\mathbf{neq}}$; if the Jacobian is banded then pd((j1) × neq + ml + mu + ij) = Jij${\mathbf{pd}}\left(\left(j-1\right)×{\mathbf{neq}}+{\mathbf{ml}}+{\mathbf{mu}}+i-j\right)={J}_{ij}$, for max (1,jmu)imin (n,j + ml)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,j-{\mathbf{mu}}\right)\le i\le \mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(n,j+{\mathbf{ml}}\right)$; (see also in nag_lapack_dgbtrf (f07bd)).
2:     user – Any MATLAB object
10:   icom(50 + neq$50+{\mathbf{neq}}$) – int64int32nag_int array
icom contains information which is usually of no interest, but is necessary for subsequent calls. However you may find the following useful:
icom(22)${\mathbf{icom}}\left(22\right)$
The order of the method to be attempted on the next step.
icom(23)${\mathbf{icom}}\left(23\right)$
The order of the method used on the last step.
icom(26)${\mathbf{icom}}\left(26\right)$
The number of steps taken so far.
icom(27)${\mathbf{icom}}\left(27\right)$
The number of calls to res so far.
icom(28)${\mathbf{icom}}\left(28\right)$
The number of evaluations of the matrix of partial derivatives needed so far.
icom(29)${\mathbf{icom}}\left(29\right)$
The total number of error test failures so far.
icom(30)${\mathbf{icom}}\left(30\right)$
The total number of convergence test failures so far.
11:   com(lcom) – double array
lcom, the dimension of the array, must satisfy the constraint lcom40 + (maxorder + 4) × neq + neq × p + q$\mathit{lcom}\ge 40+\left(\mathit{maxorder}+4\right)×{\mathbf{neq}}+{\mathbf{neq}}×p+q$ where maxorder$\mathit{maxorder}$ is the maximum order that can be used by the integration method (see maxord in nag_ode_dae_dassl_setup (d02mw)); p = neq$p={\mathbf{neq}}$ when the Jacobian is full and p = (2 × ml + mu + 1)$p=\left(2×{\mathbf{ml}}+{\mathbf{mu}}+1\right)$ when the Jacobian is banded; and, q = (neq / (ml + mu + 1)) + 1$q=\left({\mathbf{neq}}/\left({\mathbf{ml}}+{\mathbf{mu}}+1\right)\right)+1$ when the Jacobian is to be evaluated numerically and q = 0$q=0$ otherwise.
com contains information which is usually of no interest, but is necessary for subsequent calls. However you may find the following useful:
com(3)${\mathbf{com}}\left(3\right)$
The step size to be attempted on the next step.
com(4)${\mathbf{com}}\left(4\right)$
The current value of the independent variable, i.e., the farthest point integration has reached. This will be different from t only when interpolation has been performed (itask = 3${\mathbf{itask}}=3$).

### Optional Input Parameters

1:     neq – int64int32nag_int scalar
Default: The dimension of the arrays y, ydot. (An error is raised if these dimensions are not equal.)
The number of differential-algebraic equations to be solved.
Constraint: neq1${\mathbf{neq}}\ge 1$.
2:     user – Any MATLAB object
user is not used by nag_ode_dae_dassl_gen (d02ne), but is passed to res and jac. Note that for large objects it may be more efficient to use a global variable which is accessible from the m-files than to use user.

lcom iuser ruser

### Output Parameters

1:     t – double scalar
On intermediate exit: t$t$, the current value of the independent variable.
On final exit: the value of the independent variable at which the computed solution y$y$ is returned (usually at tout).
2:     y(neq) – double array
On intermediate exit: the computed solution vector, y$y$, evaluated at t = T$t=T$.
On final exit: the computed solution vector, evaluated at t$t$ (usually t = tout$t={\mathbf{tout}}$).
3:     ydot(neq) – double array
The time derivatives y ' $y\text{'}$ of the vector y$y$ at the last integration point.
4:     rtol( : $:$) – double array
Note: the dimension of the array rtol depends on the value of itol as set in nag_ode_dae_dassl_setup (d02mw); it  must be at least neq${\mathbf{neq}}$ if itol = true${\mathbf{itol}}=\mathbf{true}$ and at least 1$1$ if itol = false${\mathbf{itol}}=\mathbf{false}$.
rtol remains unchanged unless nag_ode_dae_dassl_gen (d02ne) exits with ${\mathbf{ifail}}={\mathbf{16}}$ in which case the values may have been increased to values estimated to be appropriate for continuing the integration.
5:     atol( : $:$) – double array
Note: the dimension of the array atol depends on the value of itol as set in nag_ode_dae_dassl_setup (d02mw); it  must be at least neq${\mathbf{neq}}$ if itol = true${\mathbf{itol}}=\mathbf{true}$ and at least 1$1$ if itol = false${\mathbf{itol}}=\mathbf{false}$.
atol remains unchanged unless nag_ode_dae_dassl_gen (d02ne) exits with ${\mathbf{ifail}}={\mathbf{16}}$ in which case the values may have been increased to values estimated to be appropriate for continuing the integration.
The task performed by the integrator on successful completion or an indicator that a problem occurred during integration.
itask = 2${\mathbf{itask}}=2$
The integration to tout was successfully completed (${\mathbf{t}}={\mathbf{tout}}$) by stepping exactly to tout.
itask = 3${\mathbf{itask}}=3$
The integration to tout was successfully completed (${\mathbf{t}}={\mathbf{tout}}$) by stepping past tout. y and ydot are obtained by interpolation.
itask < 0${\mathbf{itask}}<0$
Different negative values of itask returned correspond to different failure exits. ifail should always be checked in such cases and the corrective action taken where appropriate.
itask must remain unchanged between calls to nag_ode_dae_dassl_gen (d02ne).
7:     icom(50 + neq$50+{\mathbf{neq}}$) – int64int32nag_int array
8:     com(lcom) – double array
9:     user – Any MATLAB object
10:   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_ode_dae_dassl_gen (d02ne) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
 On entry, neq < 1${\mathbf{neq}}<1$
ifail = 3${\mathbf{ifail}}=3$
 On entry, ${\mathbf{tout}}={\mathbf{t}}$, or tout is too close to t to start integration, or tout is behind t in the direction of h0 (see nag_ode_dae_dassl_setup (d02mw)),
ifail = 6${\mathbf{ifail}}=6$
 On entry, rtol(i) < 0.0${\mathbf{rtol}}\left(\mathit{i}\right)<0.0$, for i = 1,2, … ,n$\mathit{i}=1,2,\dots ,n$, where n = neq$n={\mathbf{neq}}$ when itol = 1${\mathbf{itol}}=1$ and n = 1$n=1$ otherwise, or rtol(i) = atol(i) = 0.0${\mathbf{rtol}}\left(i\right)={\mathbf{atol}}\left(i\right)=0.0$ for all relevant i$i$.
ifail = 7${\mathbf{ifail}}=7$
 On entry, atol(i) < 0.0${\mathbf{atol}}\left(\mathit{i}\right)<0.0$, for i = 1,2, … ,n$\mathit{i}=1,2,\dots ,n$, where n = neq$n={\mathbf{neq}}$ when itol = 1${\mathbf{itol}}=1$ and n = 1$n=1$ otherwise.
ifail = 8${\mathbf{ifail}}=8$
 On entry, a previous call to this function returned itask < 0${\mathbf{itask}}<0$ and ifail ≠ 0${\mathbf{ifail}}\ne 0$, but no appropriate action was taken. For example, if a call returns with ${\mathbf{ifail}}={\mathbf{15}}$ (itask = -1${\mathbf{itask}}=-1$) then a call to nag_ode_dae_dassl_cont (d02mc) must be made prior to making a continuation call to nag_ode_dae_dassl_gen (d02ne).
ifail = 12${\mathbf{ifail}}=12$
Either the initialization function nag_ode_dae_dassl_setup (d02mw) has not been called prior to the first call of this function or one of the communication arrays icom, com has become corrupted.
ifail = 13${\mathbf{ifail}}=13$
 On entry, lcom < 40 + (maxorder + 4) × neq + neq × neq$\mathit{lcom}<40+\left(\mathit{maxorder}+4\right)×{\mathbf{neq}}+{\mathbf{neq}}×{\mathbf{neq}}$, and the Jacobian is full; or lcom < 40 + (maxorder + 4) × neq + neq × (2 × ml + mu + 1) + q$\mathit{lcom}<40+\left(\mathit{maxorder}+4\right)×{\mathbf{neq}}+{\mathbf{neq}}×\left(2×{\mathbf{ml}}+{\mathbf{mu}}+1\right)+q$, and the Jacobian is banded,
where maxorder$\mathit{maxorder}$ is the maximum order of the integration method to be used, and q = (neq / (ml + mu + 1)) + 1$q=\left({\mathbf{neq}}/\left({\mathbf{ml}}+{\mathbf{mu}}+1\right)\right)+1$ when the Jacobian is to be evaluated numerically and q = 0$q=0$ otherwise.
ifail = 15${\mathbf{ifail}}=15$
The maximum number of steps (500$500$) has been taken on this call and the integration has not reached tout. The integration can proceed by calling nag_ode_dae_dassl_cont (d02mc) prior to calling nag_ode_dae_dassl_gen (d02ne) again; this will reset the step counter to zero.
ifail = 16${\mathbf{ifail}}=16$
Too much accuracy was requested for the precision of the machine. On output rtol and atol were increased by an appropriate scale factor to prevent this error exit. Try running the problem again with these scaled tolerances.
ifail = 17${\mathbf{ifail}}=17$
A purely relative tolerance was selected for a given solution component, but that solution component has become zero and a pure relative error test is impossible for this component. Perhaps an absolute tolerance requirement is also necessary for this component.
ifail = 18${\mathbf{ifail}}=18$
The error test failed repeatedly using the minimum stepsize and so the integration could not proceed. If a (nonzero) minimum stepsize was specified in a call to nag_ode_dae_dassl_setup (d02mw) then either a reduction in this minimum stepsize or in the specified tolerances should be considered.
ifail = 19${\mathbf{ifail}}=19$
The corrector step to obtain the approximate solution at the next time step repeatedly failed to converge using the minimum stepsize. This may be due to an inconsistency between the system of equations to be solved and initial conditions.
ifail = 20${\mathbf{ifail}}=20$
The iteration matrix (Jacobian) has become singular. Please check the Jacobian evaluations when this is performed analytically. Also check for invalid solution values in the call to res; these should be flagged by returning ires = -1${\mathbf{ires}}=-1$. This is probably due to an error in the analytic evaluation of the Jacobian.
ifail = 21${\mathbf{ifail}}=21$
The corrector step could not converge and the error test failed repeatedly.
ifail = 22${\mathbf{ifail}}=22$
ires was set to -1$-1$ during a call to res and the problem could not be resolved.
ifail = 23${\mathbf{ifail}}=23$
ires was set to -2$-2$ during a call to res.
ifail = 24${\mathbf{ifail}}=24$
The initial ydot could not be computed. This could happen because the initial approximation to ydot was very poor or because no ydot exists that is consistent with the initial Y$Y$.
ifail = 25${\mathbf{ifail}}=25$
Repeated occurrences of input constraint violations have been detected. This could result in a potential infinite loop.

## Accuracy

The accuracy of the numerical solution may be controlled by a careful choice of the parameters rtol and atol. You are advised to use scalar error control unless the components of the solution are expected to be poorly scaled. For the type of decaying solution typical of many stiff problems, relative error control with a small absolute error threshold will be most appropriate (that is, you are advised to choose itol = 0${\mathbf{itol}}=0$ with atol(1)${\mathbf{atol}}\left(1\right)$ small but positive).

The cost of computing a solution depends critically on the size of the differential system and to a lesser extent on the degree of stiffness of the problem. For banded systems the cost is proportional to neq × (ml + mu + 1)2 ${\mathbf{neq}}×{\left({\mathbf{ml}}+{\mathbf{mu}}+1\right)}^{2}$, while for full systems the cost is proportional to neq3${{\mathbf{neq}}}^{3}$. Note however that for moderately sized problems which are only mildly nonlinear the cost may be dominated by factors proportional to neq × (ml + mu + 1)${\mathbf{neq}}×\left({\mathbf{ml}}+{\mathbf{mu}}+1\right)$ and neq2${{\mathbf{neq}}}^{2}$ respectively.

## Example

```function nag_ode_dae_dassl_gen_example
neq = 3;
maxord = 5;
mu = 2;
ml = 1;
lcom = 40+(maxord+4)*neq+(2*ml+mu+1)*neq+2*(neq/(ml+mu+1)+1)
itol = int64(1);
rtol = [1e-3; 1e-3; 1e-3];
atol = [1e-6; 1e-6; 1e-6];
ydot = zeros(neq,1);
% Set initial values
y    = [1; 0; 0];
% Initialize the problem, specifying that the Jacobian is to be
% evaluated analytically using the provided routine jac.
jceval = 'Analytic';
hmax = 0;
ho = 0;
t = 0;
tout = 0.02;
[icom, com, ifail] = ...
nag_ode_dae_dassl_setup(int64(neq), int64(maxord), jceval, hmax, ho, itol, int64(lcom));

% Specify that the Jacobian is banded
if ifail == 0
[icom, ifail] = nag_ode_dae_dassl_linalg(int64(neq), int64(ml), int64(mu), icom);
end

% Use the user parameter to pass the band dimensions through to jac.
% An alternative would be to hard code values for ml and mu in jac.
user = {ml, mu};

fprintf('\n    t            y(1)        y(2)        y(3)   \n');
fprintf(' %8.4f   %12.6f %12.6f %12.6f\n', t, y);

% Obtain the solution at 5 equally spaced values of T.
for j = 1:5
if ifail == 0
[t, y, ydot, rtol, atol, itask, icom, com, user, ifail] = ...
nag_ode_dae_dassl_gen(t, tout, y, ydot, rtol, atol, itask, @res, @jac, ...
icom, com, 'user', user);
fprintf(' %8.4f   %12.6f %12.6f %12.6f\n', t, y);
tout = tout + 0.02;
icom = nag_ode_dae_dassl_cont(icom);
end
end

function [pd, user] = jac(neq, t, y, ydot, pd, cj, user)
ml = user{1};
mu = user{2};

stride = 2*ml+mu+1;
% Main diagonal pdfull(i,i), i=1,neq
md = mu + ml + 1;
pd(md) = -0.04 - cj;
pd(md+stride) = -1.0e4*y(3) - 6.0e7*y(2) - cj;
pd(md+2*stride) = -cj;
% 1 sub-diagonal pdfull(i+1:i), i=1,neq-1
ms = md + 1;
pd(ms) = 0.04;
pd(ms+stride) = 6.0e7*y(2);
% First super-diagonal pdfull(i-1,i), i=2, neq
ms = md - 1;
pd(ms+stride) = 1.0e4*y(3);
pd(ms+2*stride) = -1.0e4*y(2);
% Second super-diagonal pdfull(i-2,i), i=3, neq
ms = md - 2;
pd(ms+2*stride) = 1.0e4*y(2);

function [r, ires, user] = res(neq, t, y, ydot, ires, user)
r = zeros(neq, 1);
r(1) = -0.04*y(1) + 1.0e4*y(2)*y(3) - ydot(1);
r(2) = 0.04*y(1) - 1.0e4*y(2)*y(3) - 3.0e7*y(2)*y(2) - ydot(2);
r(3) = 3.0e7*y(2)*y(2) - ydot(3);
```
```

lcom =

85.5000

t            y(1)        y(2)        y(3)
0.0000       1.000000     0.000000     0.000000
0.0200       0.999204     0.000036     0.000760
0.0400       0.998415     0.000036     0.001549
0.0600       0.997631     0.000036     0.002333
0.0800       0.996852     0.000036     0.003112
0.1000       0.996080     0.000036     0.003884

```
```function d02ne_example
neq = 3;
maxord = 5;
mu = 2;
ml = 1;
lcom = 40+(maxord+4)*neq+(2*ml+mu+1)*neq+2*(neq/(ml+mu+1)+1)
itol = int64(1);
rtol = [1e-3; 1e-3; 1e-3];
atol = [1e-6; 1e-6; 1e-6];
ydot = zeros(neq,1);
% Set initial values
y    = [1; 0; 0];
% Initialize the problem, specifying that the Jacobian is to be
% evaluated analytically using the provided routine jac.
jceval = 'Analytic';
hmax = 0;
ho = 0;
t = 0;
tout = 0.02;
[icom, com, ifail] = d02mw(int64(neq), int64(maxord), jceval, hmax, ho, itol, int64(lcom));

% Specify that the Jacobian is banded
if ifail == 0
[icom, ifail] = d02np(int64(neq), int64(ml), int64(mu), icom);
end

% Use the user parameter to pass the band dimensions through to jac.
% An alternative would be to hard code values for ml and mu in jac.
user = {ml, mu};

fprintf('\n    t            y(1)        y(2)        y(3)   \n');
fprintf(' %8.4f   %12.6f %12.6f %12.6f\n', t, y);

% Obtain the solution at 5 equally spaced values of T.
for j = 1:5
if ifail == 0
[t, y, ydot, rtol, atol, itask, icom, com, user, ifail] = ...
d02ne(t, tout, y, ydot, rtol, atol, itask, @res, @jac, ...
icom, com, 'user', user);
fprintf(' %8.4f   %12.6f %12.6f %12.6f\n', t, y);
tout = tout + 0.02;
icom = d02mc(icom);
end
end

function [pd, user] = jac(neq, t, y, ydot, pd, cj, user)
ml = user{1};
mu = user{2};

stride = 2*ml+mu+1;
% Main diagonal pdfull(i,i), i=1,neq
md = mu + ml + 1;
pd(md) = -0.04 - cj;
pd(md+stride) = -1.0e4*y(3) - 6.0e7*y(2) - cj;
pd(md+2*stride) = -cj;
% 1 sub-diagonal pdfull(i+1:i), i=1,neq-1
ms = md + 1;
pd(ms) = 0.04;
pd(ms+stride) = 6.0e7*y(2);
% First super-diagonal pdfull(i-1,i), i=2, neq
ms = md - 1;
pd(ms+stride) = 1.0e4*y(3);
pd(ms+2*stride) = -1.0e4*y(2);
% Second super-diagonal pdfull(i-2,i), i=3, neq
ms = md - 2;
pd(ms+2*stride) = 1.0e4*y(2);

function [r, ires, user] = res(neq, t, y, ydot, ires, user)
r = zeros(neq, 1);
r(1) = -0.04*y(1) + 1.0e4*y(2)*y(3) - ydot(1);
r(2) = 0.04*y(1) - 1.0e4*y(2)*y(3) - 3.0e7*y(2)*y(2) - ydot(2);
r(3) = 3.0e7*y(2)*y(2) - ydot(3);
```
```

lcom =

85.5000

t            y(1)        y(2)        y(3)
0.0000       1.000000     0.000000     0.000000
0.0200       0.999204     0.000036     0.000760
0.0400       0.998415     0.000036     0.001549
0.0600       0.997631     0.000036     0.002333
0.0800       0.996852     0.000036     0.003112
0.1000       0.996080     0.000036     0.003884