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_pde_2d_gen_order2_rectilinear_extractgrid (d03rz)

## Purpose

nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz) is designed to be used in conjunction with nag_pde_2d_gen_order2_rectilinear (d03rb). It can be called from the monitr to obtain the number of grid points and their (x,y)$\left(x,y\right)$ coordinates on a solution grid.

## Syntax

[npts, x, y, ifail] = d03rz(level, nlev, xmin, ymin, dxb, dyb, lgrid, istruc, lenxy)
[npts, x, y, ifail] = nag_pde_2d_gen_order2_rectilinear_extractgrid(level, nlev, xmin, ymin, dxb, dyb, lgrid, istruc, lenxy)

## Description

nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz) extracts the number of grid points and their (x,y)$\left(x,y\right)$ coordinates on a specific solution grid produced by nag_pde_2d_gen_order2_rectilinear (d03rb). It must be called only from within the monitr. The parameters nlev, xmin, ymin, dxb, dyb, lgrid and istruc to monitr must be passed unchanged to nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz).

None.

## Parameters

### Compulsory Input Parameters

1:     level – int64int32nag_int scalar
The grid level at which the coordinates are required.
Constraint: 1levelnlev$1\le {\mathbf{level}}\le {\mathbf{nlev}}$.
2:     nlev – int64int32nag_int scalar
3:     xmin – double scalar
4:     ymin – double scalar
5:     dxb – double scalar
6:     dyb – double scalar
7:     lgrid( : $:$) – int64int32nag_int array
8:     istruc( : $:$) – int64int32nag_int array
Note: the dimension of the array lgrid must be at least nlev${\mathbf{nlev}}$.
The dimension of the array istruc must be at least + 2 × nrows + npts + 1${\mathbf{lgrid}}\left({\mathbf{nlev}}\right)+2×\mathit{nrows}+{\mathbf{npts}}+1$ where nrows$\mathit{nrows}$ is stored in istruc(lgrid(nlev))${\mathbf{istruc}}\left({\mathbf{lgrid}}\left({\mathbf{nlev}}\right)\right)$ and is the number of rows in the grid at level nlev.
nlev, xmin, ymin, dxb, dyb, lgrid and istruc as supplied to monitr must be passed unchanged to nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz).
9:     lenxy – int64int32nag_int scalar
The dimension of the arrays x and y as declared in the (sub)program from which nag_pde_2d_gen_order2_rectilinear_extractgrid (d03rz) is called.
Constraint: ${\mathbf{lenxy}}\ge {\mathbf{npts}}$.

None.

None.

### Output Parameters

1:     npts – int64int32nag_int scalar
The number of grid points in the grid level level.
2:     x(lenxy) – double array
3:     y(lenxy) – double array
x(i)${\mathbf{x}}\left(\mathit{i}\right)$ and y(i)${\mathbf{y}}\left(\mathit{i}\right)$ contain the (x,y)$\left(x,y\right)$ coordinates respectively of the i$\mathit{i}$th grid point, for i = 1,2,,npts$\mathit{i}=1,2,\dots ,{\mathbf{npts}}$.
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$
 On entry, level < 1${\mathbf{level}}<1$, or ${\mathbf{level}}>{\mathbf{nlev}}$.
ifail = 2${\mathbf{ifail}}=2$
The dimension of the arrays x and y is too small for the requested grid level, i.e., ${\mathbf{lenxy}}<{\mathbf{npts}}$.

Not applicable.

None.

## Example

```function nag_pde_2d_gen_order2_rectilinear_extractgrid_example
ts = 0;
twant = [0.25; 1];
dt = [0.001;
1e-07;
0];
tols = 0.1;
tolt = 0.05;
opti = [int64(5);0;0;0];
optr = [1, 1;
1, 1;
1, 1];
rwk = zeros(426000, 1);
iwk = zeros(117037, 1, 'int64');
lenlwk = int64(6000);
itrace = int64(0);
ind = int64(0);
global iout;
for iout = 1:2
tout = twant(iout);
[ts, dt, rwk, iwk, ind, ifail] = ...
nag_pde_2d_gen_order2_rectilinear(ts, tout, dt, tols, tolt, @inidom, @pdedef, ...
@bndary, @pdeiv, @monitr, opti, optr, ...
rwk, iwk, lenlwk, itrace, ind);
fprintf('\nStatistics\n');
fprintf('Time = %8.4f\n', ts);
fprintf('Total number of accepted timesteps = %d\n', iwk(1));
fprintf('Total number of rejected timesteps = %d\n', iwk(2));
fprintf('\n             T o t a l  n u m b e r  o f    \n');
fprintf('          Residual  Jacobian    Newton  Lin sys\n');
fprintf('             evals     evals     iters    iters\n');
fprintf(' At level \n');
maxlev = double(opti(1));
for j = 1:maxlev
if (iwk(j+2) ~= 0)
fprintf('%6d %10d %10d %10d %10d\n', j, iwk(j+2), iwk(j+2+maxlev), ...
iwk(j+2+2*maxlev), iwk(j+2+3*maxlev) );
end
end
fprintf('\n             M a x i m u m  n u m b e r  o f\n');
fprintf('              Newton iters    Lin sys iters \n');
fprintf(' At level \n');
for j = 1:maxlev
if (iwk(j+2) ~= 0)
fprintf('%6d%14d%14d\n', j, iwk(j+2+4*maxlev), iwk(j+2+5*maxlev) );
end
end

end

function [res] = ...
bndary(npts, npde, t, x, y, u, ut, ux, uy, nbnds, nbpts, llbnd, ilbnd, lbnd, res)

epsilon = 1e-3;

for k = double(llbnd(1)):double(nbpts)
i = lbnd(k);
a = (-4.0d0*x(i)+4.0d0*y(i)-t)/(32.0d0*epsilon);
if (a <= 0.0d0)
res(i,1) = u(i,1) - (0.75d0-0.25d0/(1.0d0+exp(a)));
res(i,2) = u(i,2) - (0.75d0+0.25d0/(1.0d0+exp(a)));
else
res(i,1) = u(i,1) - (0.75d0-0.25d0*exp(-a)/(exp(-a)+1.0d0));
res(i,2) = u(i,2) - (0.75d0+0.25d0*exp(-a)/(exp(-a)+1.0d0));
end
end

function [xmin, xmax, ymin, ymax, nx, ny, npts, nrows, nbnds, nbpts, ...
lrow, irow, icol, llbnd, ilbnd, lbnd, ierr] = inidom(maxpts, ierr)

nrows = int64(11);
npts = int64(105);
nbnds = int64(28);
nbpts = int64(72);

lrow = zeros(nrows, 1, 'int64');
irow = zeros(nrows, 1, 'int64');
icol = zeros(npts, 1, 'int64');
llbnd = zeros(nbnds, 1, 'int64');
ilbnd = zeros(nbnds, 1, 'int64');
lbnd = zeros(nbpts, 1, 'int64');

icold = [int64(0),1,2,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,9, ...
10,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,8,9,10,0,1,2,3,4,5, ...
6,7,8,9,10,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,9,10, ...
0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8];

ilbndd = [int64(1),2,3,4,1,4,1,2,3,4,3,4,1,2,12,23,34,41,14,41,12, ...
23,34,41,43,14,21,32];

irowd = [int64(0),1,2,3,4,5,6,7,8,9,10];

lbndd = [int64(2),4,15,26,37,46,57,68,79,88,98,99,100,101,102,103, ...
104,96,86,85,84,83,82,70,59,48,39,28,17,6,8,9,10,11,12,13, ...
18,29,40,49,60,72,73,74,75,76,77,67,56,45,36,25,33,32,42,52, ...
53,43,1,97,105,87,81,3,7,71,78,14,31,51,54,34];

llbndd = [int64(1),2,11,18,19,24,31,37,42,48,53,55,56,58,59,60,61,62, ...
63,64,65,66,67,68,69,70,71,72];

lrowd = [int64(1),4,15,26,37,46,57,68,79,88,97];

nx = int64(11);
ny = int64(11);

% check maxpts against rough estimate of npts
if (maxpts < double(nx)*double(ny))
ierr = -1;
return;
end

xmin = 0.0d0;
ymin = 0.0d0;
xmax = 1.0d0;
ymax = 1.0d0;

for i = 1:double(nrows)
lrow(i) = lrowd(i);
irow(i) = irowd(i);
end

for i = 1:double(nbnds)
llbnd(i) = llbndd(i);
ilbnd(i) = ilbndd(i);
end

for i = 1:double(nbpts)
lbnd(i) = lbndd(i);
end

for i = 1:double(npts)
icol(i) = icold(i);
end

function [res] = pdedef(npts, npde, t, x, y, u, ut, ux, uy, uxx, uxy, uyy)
res = zeros(npts, npde);

epsilon = 1e-3;
for i = 1:double(npts)
res(i,1) = ut(i,1)-(-u(i,1)*ux(i,1)-u(i,2)*uy(i,1)+epsilon*(uxx(i,1)+uyy(i,1)));
res(i,2) = ut(i,2)-(-u(i,1)*ux(i,2)-u(i,2)*uy(i,2)+epsilon*(uxx(i,2)+uyy(i,2)));
end

function [u] = pdeiv(npts, npde, t, x, y)
u = zeros(npts, npde);

epsilon = 1e-3;

for i = 1:double(npts)
a = (-4.0d0*x(i)+4.0d0*y(i)-t)/(32.0d0*epsilon);
if (a <= 0.0d0)
u(i,1) = 0.75d0 - 0.25d0/(1.0d0+exp(a));
u(i,2) = 0.75d0 + 0.25d0/(1.0d0+exp(a));
else
u(i,1) = 0.75d0 - 0.25d0*exp(-a)/(exp(-a)+1.0d0);
u(i,2) = 0.75d0 + 0.25d0*exp(-a)/(exp(-a)+1.0d0);
end
end

function [ierr] = ...
monitr(npde, t, dt, dtnew, tlast, nlev, xmin, ymin, dxb, dyb, lgrid, ...
istruc, lsol, sol, ierr)

lenxy=int64(2500);

global iout;

if (tlast && iout == 2)
for level = 1:double(nlev)
ipsol = lsol(level);

% Get grid information
[npts, x, y, ifail] = ...
nag_pde_2d_gen_order2_rectilinear_extractgrid(int64(level), nlev,  ...
xmin, ymin, dxb, dyb, lgrid, istruc, lenxy);
if (ifail ~= 0)
ierr = 1;
elseif (level == 1)
% Get exact solution
[uex] = pdeiv(npts, npde, t, x, y);

fprintf('\n Solution at every 2nd grid point in level 1 at time %8.4f\n', t);
fprintf('       x          y        approx u     exact u    approx v     exact v\n');
for ipt = 1:2:double(npts)
fprintf('%11.4E %11.4e %11.4e %11.4e %11.4e %11.4e\n', x(ipt), ...
y(ipt), sol(double(ipsol)+double(ipt)), uex(ipt,1), ...
sol(double(ipsol)+double(npts)+double(ipt)), uex(ipt,2));
end
fprintf('\n');
end
end
end
```
```

Statistics
Time =   0.2500
Total number of accepted timesteps = 14
Total number of rejected timesteps = 0

T o t a l  n u m b e r  o f
Residual  Jacobian    Newton  Lin sys
evals     evals     iters    iters
At level
1        196         14         28         14
2        196         14         28         22
3        196         14         28         25
4        196         14         28         31
5        141         10         21         29

M a x i m u m  n u m b e r  o f
Newton iters    Lin sys iters
At level
1             2             1
2             2             1
3             2             1
4             2             2
5             3             2

Solution at every 2nd grid point in level 1 at time   1.0000
x          y        approx u     exact u    approx v     exact v
0.0000E+00  0.0000e+00  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
2.0000E-01  0.0000e+00  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
1.0000E-01  1.0000e-01  5.0022e-01  5.0000e-01  9.9978e-01  1.0000e+00
3.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
5.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
7.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
9.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
0.0000E+00  2.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
2.0000E-01  2.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
4.0000E-01  2.0000e-01  5.0011e-01  5.0000e-01  9.9989e-01  1.0000e+00
6.0000E-01  2.0000e-01  4.9994e-01  5.0000e-01  1.0001e+00  1.0000e+00
8.0000E-01  2.0000e-01  4.9999e-01  5.0000e-01  1.0000e+00  1.0000e+00
1.0000E+00  2.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
1.0000E-01  3.0000e-01  4.9997e-01  5.0048e-01  1.0000e+00  9.9952e-01
3.0000E-01  3.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
5.0000E-01  3.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
7.0000E-01  3.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
9.0000E-01  3.0000e-01  5.0003e-01  5.0000e-01  9.9997e-01  1.0000e+00
0.0000E+00  4.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
2.0000E-01  4.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
4.0000E-01  4.0000e-01  5.0018e-01  5.0000e-01  9.9982e-01  1.0000e+00
8.0000E-01  4.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
1.0000E+00  4.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
1.0000E-01  5.0000e-01  7.4998e-01  7.5000e-01  7.5002e-01  7.5000e-01
3.0000E-01  5.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
5.0000E-01  5.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
7.0000E-01  5.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
9.0000E-01  5.0000e-01  5.0005e-01  5.0000e-01  9.9995e-01  1.0000e+00
0.0000E+00  6.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
2.0000E-01  6.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
4.0000E-01  6.0000e-01  4.9998e-01  5.0048e-01  1.0000e+00  9.9952e-01
6.0000E-01  6.0000e-01  4.9991e-01  5.0000e-01  1.0001e+00  1.0000e+00
8.0000E-01  6.0000e-01  4.9983e-01  5.0000e-01  1.0002e+00  1.0000e+00
1.0000E+00  6.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
1.0000E-01  7.0000e-01  7.5002e-01  7.5000e-01  7.4998e-01  7.5000e-01
3.0000E-01  7.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
5.0000E-01  7.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
7.0000E-01  7.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
9.0000E-01  7.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
0.0000E+00  8.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
2.0000E-01  8.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
4.0000E-01  8.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
6.0000E-01  8.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
8.0000E-01  8.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
1.0000E-01  9.0000e-01  7.5001e-01  7.5000e-01  7.4999e-01  7.5000e-01
3.0000E-01  9.0000e-01  7.5002e-01  7.5000e-01  7.4998e-01  7.5000e-01
5.0000E-01  9.0000e-01  7.5002e-01  7.5000e-01  7.4998e-01  7.5000e-01
7.0000E-01  9.0000e-01  4.9994e-01  5.0048e-01  1.0001e+00  9.9952e-01
0.0000E+00  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
2.0000E-01  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
4.0000E-01  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
6.0000E-01  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
8.0000E-01  1.0000e+00  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01

Statistics
Time =   1.0000
Total number of accepted timesteps = 45
Total number of rejected timesteps = 0

T o t a l  n u m b e r  o f
Residual  Jacobian    Newton  Lin sys
evals     evals     iters    iters
At level
1        630         45         90         45
2        630         45         90         78
3        630         45         90         87
4        630         45         90        124
5        575         41         83        122

M a x i m u m  n u m b e r  o f
Newton iters    Lin sys iters
At level
1             2             1
2             2             1
3             2             1
4             2             2
5             3             2

```
```function d03rz_example
ts = 0;
twant = [0.25; 1];
dt = [0.001;
1e-07;
0];
tols = 0.1;
tolt = 0.05;
opti = [int64(5);0;0;0];
optr = [1, 1;
1, 1;
1, 1];
rwk = zeros(426000, 1);
iwk = zeros(117037, 1, 'int64');
lenlwk = int64(6000);
itrace = int64(0);
ind = int64(0);
global iout;
for iout = 1:2
tout = twant(iout);
[ts, dt, rwk, iwk, ind, ifail] = ...
d03rb(ts, tout, dt, tols, tolt, @inidom, @pdedef, ...
@bndary, @pdeiv, @monitr, opti, optr, ...
rwk, iwk, lenlwk, itrace, ind);
fprintf('\nStatistics\n');
fprintf('Time = %8.4f\n', ts);
fprintf('Total number of accepted timesteps = %d\n', iwk(1));
fprintf('Total number of rejected timesteps = %d\n', iwk(2));
fprintf('\n             T o t a l  n u m b e r  o f    \n');
fprintf('          Residual  Jacobian    Newton  Lin sys\n');
fprintf('             evals     evals     iters    iters\n');
fprintf(' At level \n');
maxlev = double(opti(1));
for j = 1:maxlev
if (iwk(j+2) ~= 0)
fprintf('%6d %10d %10d %10d %10d\n', j, iwk(j+2), iwk(j+2+maxlev), ...
iwk(j+2+2*maxlev), iwk(j+2+3*maxlev) );
end
end
fprintf('\n             M a x i m u m  n u m b e r  o f\n');
fprintf('              Newton iters    Lin sys iters \n');
fprintf(' At level \n');
for j = 1:maxlev
if (iwk(j+2) ~= 0)
fprintf('%6d%14d%14d\n', j, iwk(j+2+4*maxlev), iwk(j+2+5*maxlev) );
end
end

end

function [res] = ...
bndary(npts, npde, t, x, y, u, ut, ux, uy, nbnds, nbpts, llbnd, ilbnd, lbnd, res)

epsilon = 1e-3;

for k = double(llbnd(1)):double(nbpts)
i = lbnd(k);
a = (-4.0d0*x(i)+4.0d0*y(i)-t)/(32.0d0*epsilon);
if (a <= 0.0d0)
res(i,1) = u(i,1) - (0.75d0-0.25d0/(1.0d0+exp(a)));
res(i,2) = u(i,2) - (0.75d0+0.25d0/(1.0d0+exp(a)));
else
res(i,1) = u(i,1) - (0.75d0-0.25d0*exp(-a)/(exp(-a)+1.0d0));
res(i,2) = u(i,2) - (0.75d0+0.25d0*exp(-a)/(exp(-a)+1.0d0));
end
end

function [xmin, xmax, ymin, ymax, nx, ny, npts, nrows, nbnds, nbpts, ...
lrow, irow, icol, llbnd, ilbnd, lbnd, ierr] = inidom(maxpts, ierr)

nrows = int64(11);
npts = int64(105);
nbnds = int64(28);
nbpts = int64(72);

lrow = zeros(nrows, 1, 'int64');
irow = zeros(nrows, 1, 'int64');
icol = zeros(npts, 1, 'int64');
llbnd = zeros(nbnds, 1, 'int64');
ilbnd = zeros(nbnds, 1, 'int64');
lbnd = zeros(nbpts, 1, 'int64');

icold = [int64(0),1,2,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,9, ...
10,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,8,9,10,0,1,2,3,4,5, ...
6,7,8,9,10,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,9,10, ...
0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8];

ilbndd = [int64(1),2,3,4,1,4,1,2,3,4,3,4,1,2,12,23,34,41,14,41,12, ...
23,34,41,43,14,21,32];

irowd = [int64(0),1,2,3,4,5,6,7,8,9,10];

lbndd = [int64(2),4,15,26,37,46,57,68,79,88,98,99,100,101,102,103, ...
104,96,86,85,84,83,82,70,59,48,39,28,17,6,8,9,10,11,12,13, ...
18,29,40,49,60,72,73,74,75,76,77,67,56,45,36,25,33,32,42,52, ...
53,43,1,97,105,87,81,3,7,71,78,14,31,51,54,34];

llbndd = [int64(1),2,11,18,19,24,31,37,42,48,53,55,56,58,59,60,61,62, ...
63,64,65,66,67,68,69,70,71,72];

lrowd = [int64(1),4,15,26,37,46,57,68,79,88,97];

nx = int64(11);
ny = int64(11);

% check maxpts against rough estimate of npts
if (maxpts < double(nx)*double(ny))
ierr = -1;
return;
end

xmin = 0.0d0;
ymin = 0.0d0;
xmax = 1.0d0;
ymax = 1.0d0;

for i = 1:double(nrows)
lrow(i) = lrowd(i);
irow(i) = irowd(i);
end

for i = 1:double(nbnds)
llbnd(i) = llbndd(i);
ilbnd(i) = ilbndd(i);
end

for i = 1:double(nbpts)
lbnd(i) = lbndd(i);
end

for i = 1:double(npts)
icol(i) = icold(i);
end

function [res] = pdedef(npts, npde, t, x, y, u, ut, ux, uy, uxx, uxy, uyy)
res = zeros(npts, npde);

epsilon = 1e-3;
for i = 1:double(npts)
res(i,1) = ut(i,1)-(-u(i,1)*ux(i,1)-u(i,2)*uy(i,1)+epsilon*(uxx(i,1)+uyy(i,1)));
res(i,2) = ut(i,2)-(-u(i,1)*ux(i,2)-u(i,2)*uy(i,2)+epsilon*(uxx(i,2)+uyy(i,2)));
end

function [u] = pdeiv(npts, npde, t, x, y)
u = zeros(npts, npde);

epsilon = 1e-3;

for i = 1:double(npts)
a = (-4.0d0*x(i)+4.0d0*y(i)-t)/(32.0d0*epsilon);
if (a <= 0.0d0)
u(i,1) = 0.75d0 - 0.25d0/(1.0d0+exp(a));
u(i,2) = 0.75d0 + 0.25d0/(1.0d0+exp(a));
else
u(i,1) = 0.75d0 - 0.25d0*exp(-a)/(exp(-a)+1.0d0);
u(i,2) = 0.75d0 + 0.25d0*exp(-a)/(exp(-a)+1.0d0);
end
end

function [ierr] = ...
monitr(npde, t, dt, dtnew, tlast, nlev, xmin, ymin, dxb, dyb, lgrid, ...
istruc, lsol, sol, ierr)

lenxy=int64(2500);

global iout;

if (tlast && iout == 2)
for level = 1:double(nlev)
ipsol = lsol(level);

% Get grid information
[npts, x, y, ifail] = ...
d03rz(int64(level), nlev, xmin, ymin, dxb, dyb, lgrid, istruc, lenxy);
if (ifail ~= 0)
ierr = 1;
elseif (level == 1)
% Get exact solution
[uex] = pdeiv(npts, npde, t, x, y);

fprintf('\n Solution at every 2nd grid point in level 1 at time %8.4f\n', t);
fprintf('       x          y        approx u     exact u    approx v     exact v\n');
for ipt = 1:2:double(npts)
fprintf('%11.4E %11.4e %11.4e %11.4e %11.4e %11.4e\n', x(ipt), ...
y(ipt), sol(double(ipsol)+double(ipt)), uex(ipt,1), ...
sol(double(ipsol)+double(npts)+double(ipt)), uex(ipt,2));
end
fprintf('\n');
end
end
end
```
```

Statistics
Time =   0.2500
Total number of accepted timesteps = 14
Total number of rejected timesteps = 0

T o t a l  n u m b e r  o f
Residual  Jacobian    Newton  Lin sys
evals     evals     iters    iters
At level
1        196         14         28         14
2        196         14         28         22
3        196         14         28         25
4        196         14         28         31
5        141         10         21         29

M a x i m u m  n u m b e r  o f
Newton iters    Lin sys iters
At level
1             2             1
2             2             1
3             2             1
4             2             2
5             3             2

Solution at every 2nd grid point in level 1 at time   1.0000
x          y        approx u     exact u    approx v     exact v
0.0000E+00  0.0000e+00  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
2.0000E-01  0.0000e+00  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
1.0000E-01  1.0000e-01  5.0022e-01  5.0000e-01  9.9978e-01  1.0000e+00
3.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
5.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
7.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
9.0000E-01  1.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
0.0000E+00  2.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
2.0000E-01  2.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
4.0000E-01  2.0000e-01  5.0011e-01  5.0000e-01  9.9989e-01  1.0000e+00
6.0000E-01  2.0000e-01  4.9994e-01  5.0000e-01  1.0001e+00  1.0000e+00
8.0000E-01  2.0000e-01  4.9999e-01  5.0000e-01  1.0000e+00  1.0000e+00
1.0000E+00  2.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
1.0000E-01  3.0000e-01  4.9997e-01  5.0048e-01  1.0000e+00  9.9952e-01
3.0000E-01  3.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
5.0000E-01  3.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
7.0000E-01  3.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
9.0000E-01  3.0000e-01  5.0003e-01  5.0000e-01  9.9997e-01  1.0000e+00
0.0000E+00  4.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
2.0000E-01  4.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
4.0000E-01  4.0000e-01  5.0018e-01  5.0000e-01  9.9982e-01  1.0000e+00
8.0000E-01  4.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
1.0000E+00  4.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
1.0000E-01  5.0000e-01  7.4998e-01  7.5000e-01  7.5002e-01  7.5000e-01
3.0000E-01  5.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
5.0000E-01  5.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
7.0000E-01  5.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
9.0000E-01  5.0000e-01  5.0005e-01  5.0000e-01  9.9995e-01  1.0000e+00
0.0000E+00  6.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
2.0000E-01  6.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
4.0000E-01  6.0000e-01  4.9998e-01  5.0048e-01  1.0000e+00  9.9952e-01
6.0000E-01  6.0000e-01  4.9991e-01  5.0000e-01  1.0001e+00  1.0000e+00
8.0000E-01  6.0000e-01  4.9983e-01  5.0000e-01  1.0002e+00  1.0000e+00
1.0000E+00  6.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
1.0000E-01  7.0000e-01  7.5002e-01  7.5000e-01  7.4998e-01  7.5000e-01
3.0000E-01  7.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
5.0000E-01  7.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
7.0000E-01  7.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
9.0000E-01  7.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
0.0000E+00  8.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
2.0000E-01  8.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
4.0000E-01  8.0000e-01  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
6.0000E-01  8.0000e-01  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01
8.0000E-01  8.0000e-01  5.0000e-01  5.0000e-01  1.0000e+00  1.0000e+00
1.0000E-01  9.0000e-01  7.5001e-01  7.5000e-01  7.4999e-01  7.5000e-01
3.0000E-01  9.0000e-01  7.5002e-01  7.5000e-01  7.4998e-01  7.5000e-01
5.0000E-01  9.0000e-01  7.5002e-01  7.5000e-01  7.4998e-01  7.5000e-01
7.0000E-01  9.0000e-01  4.9994e-01  5.0048e-01  1.0001e+00  9.9952e-01
0.0000E+00  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
2.0000E-01  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
4.0000E-01  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
6.0000E-01  1.0000e+00  7.5000e-01  7.5000e-01  7.5000e-01  7.5000e-01
8.0000E-01  1.0000e+00  5.0048e-01  5.0048e-01  9.9952e-01  9.9952e-01

Statistics
Time =   1.0000
Total number of accepted timesteps = 45
Total number of rejected timesteps = 0

T o t a l  n u m b e r  o f
Residual  Jacobian    Newton  Lin sys
evals     evals     iters    iters
At level
1        630         45         90         45
2        630         45         90         78
3        630         45         90         87
4        630         45         90        124
5        575         41         83        122

M a x i m u m  n u m b e r  o f
Newton iters    Lin sys iters
At level
1             2             1
2             2             1
3             2             1
4             2             2
5             3             2

```