BVPFD

Solves a (parameterized) system of differential equations with boundary conditions at two points, using a variable order, variable step size finite difference method with deferred corrections.

Required Arguments

FCNEQN — User-supplied SUBROUTINE to evaluate derivatives. The usage is CALL    FCNEQN (N, T, Y, P, DYDT), where
                                N – Number of differential equations.   (Input)
                                T – Independent variable, t.   (Input)
                                Y – Array of size N containing the dependent variable values, y(t).
                                (Input)
                                P – Continuation parameter, p.   (Input)
                                See Comment 3.
                                DYDT – Array of size N containing the derivatives yʹ (t).   (Output)
           The name FCNEQN must be declared EXTERNAL in the calling program.

FCNJAC — User-supplied SUBROUTINE to evaluate the Jacobian. The usage is CALL    FCNJAC (N, T, Y, P, DYPDY), where
                                N – Number of differential equations.   (Input)
                                T – Independent variable, t.   (Input)
                                Y – Array of size N containing the dependent variable values.   (Input)
                                P – Continuation parameter, p.   (Input)
                                See Comments 3.
                                DYPDY – N by N array containing the partial derivatives ai, j = ∂ fi ∕ ∂ yj
                                evaluated at (t, y). The values ai,j are returned in DYPDY(i, j).
                                (Output)                                                                                                                                                     The name FCNJAC must be declared EXTERNAL in the calling program.

FCNBC — User-supplied SUBROUTINE to evaluate the boundary conditions. The usage is               CALL FCNBC (N, YLEFT, YRIGHT, P, H), where
                                N – Number of differential equations.   (Input)
                                YLEFT – Array of size N containing the values of the dependent
                                variable at the left endpoint.   (Input)
                                YRIGHT – Array of size N containing the values of the dependent
                                variable at the right endpoint.   (Input)
                                P – Continuation parameter, p.   (Input)
                                See Comment 3.
                                H – Array of size N containing the boundary condition residuals.
                                (Output)

The boundary conditions are defined by hi = 0; for i = 1, …, N. The left endpoint conditions must be defined first, then, the conditions involving both endpoints, and finally the right endpoint conditions.

The name FCNBC must be declared EXTERNAL in the calling program.

FCNPEQ — User-supplied SUBROUTINE to evaluate the derivative of yʹ with respect to the parameter p. The usage is

                                            CALL FCNPEQ (N, T, Y, P, DYPDP), where
                                N – Number of differential equations.   (Input)
                                T – Dependent variable, t.   (Input)
                                Y – Array of size N containing the dependent variable values.   (Input)
                                P – Continuation parameter, p.   (Input)
                                See Comment 3.
                                DYPDP – Array of size N containing the derivative of yʹ
                                evaluated at (t, y).
                                (Output)

The name FCNPEQ must be declared EXTERNAL in the calling program.

FCNPBC — User-supplied SUBROUTINE to evaluate the derivative of the boundary           conditions with respect to the parameter p. The usage is
           CALL FCNPBC (N, YLEFT, YRIGHT, P, H), where
                                N – Number of differential equations.   (Input)
                                YLEFT – Array of size N containing the values of the dependent
                                variable at the left endpoint.   (Input)
                                YRIGHT – Array of size N containing the values of the dependent
                                variable at the right endpoint.   (Input)
                                P – Continuation parameter, p.   (Input)
                                See Comment 3.
                                H – Array of size N containing the derivative of fi with respect to p.
                                (Output)

The name FCNPBC must be declared EXTERNAL in the calling program.

NLEFT — Number of initial conditions. (Input)
The value NLEFT must be greater than or equal to zero and less than N.

NCUPBC — Number of coupled boundary conditions. (Input)
The value NLEFT + NCUPBC must be greater than zero and less than or equal to N.

TLEFT — The left endpoint. (Input)

TRIGHT — The right endpoint. (Input)

PISTEP — Initial increment size for p. (Input)
If this value is zero, continuation will not be used in this problem. The routines FCNPEQ and FCNPBC will not be called.

TOL — Relative error control parameter. (Input)
The computations stop when ABS(ERROR(J, I))/MAX(ABS(Y(J, I)), 1.0).LT.TOL for all J = 1, …, N and I = 1, …, NGRID. Here, ERROR(J, I) is the estimated error in Y(J, I).

TINIT — Array of size NINIT containing the initial grid points. (Input)

YINIT — Array of size N by NINIT containing an initial guess for the values of Y at the points in TINIT. (Input)

LINEAR — Logical .TRUE. if the differential equations and the boundary conditions are linear. (Input)

MXGRID — Maximum number of grid points allowed. (Input)

NFINAL — Number of final grid points, including the endpoints. (Output)

TFINAL — Array of size MXGRID containing the final grid points. (Output)
Only the first NFINAL points are significant.

YFINAL — Array of size N by MXGRID containing the values of Y at the points in TFINAL. (Output)

ERREST — Array of size N. (Output)
ERREST(J) is the estimated error in Y(J).

Optional Arguments

N — Number of differential equations. (Input)
Default: N = size (YINIT,1).

NINIT — Number of initial grid points, including the endpoints. (Input)
It must be at least 4.
Default: NINIT = size (TINIT,1).

LDYINI — Leading dimension of YINIT exactly as specified in the dimension statement of the calling program. (Input)
Default: LDYINI = size (YINIT,1).

PRINT — Logical .TRUE. if intermediate output is to be printed. (Input)
Default: PRINT = .FALSE.

LDYFIN — Leading dimension of YFINAL exactly as specified in the dimension statement of the calling program. (Input)
Default: LDYFIN = size (YFINAL,1).

FORTRAN 90 Interface

Generic: CALL BVPFD (FCNEQN, FCNJAC, FCNBC, FCNPEQ, FCNPBC, NLEFT, NCUPBC, TLEFT, TRIGHT, PISTEP, TOL, TINIT, YINIT, LINEAR, MXGRID, NFINAL, TFINAL, YFINAL, ERREST [,…])

Specific: The specific interface names are S_BVPFD and D_BVPFD.

FORTRAN 77 Interface

Single: CALL BVPFD (FCNEQN, FCNJAC, FCNBC, FCNPEQ, FCNPBC, N, NLEFT, NCUPBC, TLEFT, TRIGHT, PISTEP, TOL, NINIT, TINIT, YINIT, LDYINI, LINEAR, PRINT, MXGRID, NFINAL, TFINAL, YFINAL, LDYFIN, ERREST)

Double: The double precision name is DBVPFD.

Description

The routine BVPFD is based on the subprogram PASVA3 by M. Lentini and V. Pereyra (see Pereyra 1978). The basic discretization is the trapezoidal rule over a nonuniform mesh. This mesh is chosen adaptively, to make the local error approximately the same size everywhere. Higher-order discretizations are obtained by deferred corrections. Global error estimates are produced to control the computation. The resulting nonlinear algebraic system is solved by Newton's method with step control. The linearized system of equations is solved by a special form of Gauss elimination that preserves the sparseness.

Comments

1. Workspace may be explicitly provided, if desired, by use of B2PFD/DB2PFD. The reference is:

CALL B2PFD (FCNEQN, FCNJAC, FCNBC, FCNPEQ, FCNPBC, N, NLEFT, NCUPBC, TLEFT, TRIGHT, PISTEP, TOL, NINIT, TINIT, YINIT, LDYINI, LINEAR, PRINT, MXGRID, NFINAL, TFINAL, YFINAL, LDYFIN, ERREST, RWORK, IWORK)

The additional arguments are as follows:

RWORK — Floating-point work array of size N(3N * MXGRID + 4N + 1) + MXGRID * (7N + 2).

IWORK — Integer work array of size 2N * MXGRID + N + MXGRID.

2. Informational errors

Type Code

4 1 More than MXGRID grid points are needed to solve the problem.

4 2 Newton's method diverged.

3 3 Newton's method reached roundoff error level.

3. If the value of PISTEP is greater than zero, then the routine BVPFD assumes that the user has embedded the problem into a one-parameter family of problems:

yʹ = yʹ(t, y, p)

h(ytleft, ytright, p) = 0

such that for p = 0 the problem is simple. For p = 1, the original problem is recovered. The routine BVPFD automatically attempts to increment from p = 0 to p = 1. The value PISTEP is the beginning increment used in this continuation. The increment will usually be changed by routine BVPFD, but an arbitrary minimum of 0.01 is imposed.

4. The vectors TINIT and TFINAL may be the same.

5. The arrays YINIT and YFINAL may be the same.

Example 1

This example solves the third-order linear equation

subject to the boundary conditions y(0) = y(2π) and yʹ(0) = yʹ(2π) = 1. (Its solution is y = sin t.) To use BVPFD, the problem is reduced to a system of first-order equations by defining
y1 = y, y2= yʹ and y3 = y″. The resulting system is

Note that there is one boundary condition at the left endpoint t = 0 and one boundary condition coupling the left and right endpoints. The final boundary condition is at the right endpoint. The total number of boundary conditions must be the same as the number of equations (in this case 3).

Note that since the parameter p is not used in the call to BVPFD, the routines FCNPEQ and FCNPBC are not needed. Therefore, in the call to BVPFD, FCNEQN and FCNBC were used in place of FCNPEQ and FCNPBC.

USE BVPFD_INT

USE UMACH_INT

USE CONST_INT

IMPLICIT NONE

! SPECIFICATIONS FOR PARAMETERS

INTEGER LDYFIN, LDYINI, MXGRID, NEQNS, NINIT

PARAMETER (MXGRID=45, NEQNS=3, NINIT=10, LDYFIN=NEQNS, &

LDYINI=NEQNS)

! SPECIFICATIONS FOR LOCAL VARIABLES

INTEGER I, J, NCUPBC, NFINAL, NLEFT, NOUT

REAL ERREST(NEQNS), PISTEP, TFINAL(MXGRID), TINIT(NINIT), &

TLEFT, TOL, TRIGHT, YFINAL(LDYFIN,MXGRID), &

YINIT(LDYINI,NINIT)

LOGICAL LINEAR, PRINT

! SPECIFICATIONS FOR INTRINSICS

INTRINSIC FLOAT

REAL FLOAT

! SPECIFICATIONS FOR SUBROUTINES

! SPECIFICATIONS FOR FUNCTIONS

EXTERNAL FCNBC, FCNEQN, FCNJAC

! Set parameters

NLEFT = 1

NCUPBC = 1

TOL = .001

TLEFT = 0.0

TRIGHT = CONST('PI')

TRIGHT = 2.0*TRIGHT

PISTEP = 0.0

PRINT = .FALSE.

LINEAR = .TRUE.

! Define TINIT

DO 10 I=1, NINIT

TINIT(I) = TLEFT + (I-1)*(TRIGHT-TLEFT)/FLOAT(NINIT-1)

10 CONTINUE

! Set YINIT to zero

YINIT = 0.0E0

! Solve problem

CALL BVPFD (FCNEQN, FCNJAC, FCNBC, FCNEQN, FCNBC, NLEFT, &

NCUPBC, TLEFT, TRIGHT, PISTEP, TOL, TINIT, &

YINIT, LINEAR, MXGRID, NFINAL, &

TFINAL, YFINAL, ERREST)

! Print results

CALL UMACH (2, NOUT)

WRITE (NOUT,99997)

WRITE (NOUT,99998) (I,TFINAL(I),(YFINAL(J,I),J=1,NEQNS),I=1, &

NFINAL)

WRITE (NOUT,99999) (ERREST(J),J=1,NEQNS)

99997 FORMAT (4X, 'I', 7X, 'T', 14X, 'Y1', 13X, 'Y2', 13X, 'Y3')

99998 FORMAT (I5, 1P4E15.6)

99999 FORMAT (' Error estimates', 4X, 1P3E15.6)

END

SUBROUTINE FCNEQN (NEQNS, T, Y, P, DYDX)

! SPECIFICATIONS FOR ARGUMENTS

INTEGER NEQNS

REAL T, P, Y(NEQNS), DYDX(NEQNS)

! SPECIFICATIONS FOR INTRINSICS

INTRINSIC SIN

REAL SIN

! Define PDE

DYDX(1) = Y(2)

DYDX(2) = Y(3)

DYDX(3) = 2.0*Y(3) - Y(2) + Y(1) + SIN(T)

RETURN

END

SUBROUTINE FCNJAC (NEQNS, T, Y, P, DYPDY)

! SPECIFICATIONS FOR ARGUMENTS

INTEGER NEQNS

REAL T, P, Y(NEQNS), DYPDY(NEQNS,NEQNS)

! Define d(DYDX)/dY

DYPDY(1,1) = 0.0

DYPDY(1,2) = 1.0

DYPDY(1,3) = 0.0

DYPDY(2,1) = 0.0

DYPDY(2,2) = 0.0

DYPDY(2,3) = 1.0

DYPDY(3,1) = 1.0

DYPDY(3,2) = -1.0

DYPDY(3,3) = 2.0

RETURN

END

SUBROUTINE FCNBC (NEQNS, YLEFT, YRIGHT, P, F)

! SPECIFICATIONS FOR ARGUMENTS

INTEGER NEQNS

REAL P, YLEFT(NEQNS), YRIGHT(NEQNS), F(NEQNS)

! Define boundary conditions

F(1) = YLEFT(2) - 1.0

F(2) = YLEFT(1) - YRIGHT(1)

F(3) = YRIGHT(2) - 1.0

RETURN

END

Output

I       T              Y1             Y2             Y3
1   0.000000E+00 -1.123191E-04   1.000000E+00   6.242319E-05
2   3.490659E-01   3.419107E-01   9.397087E-01 -3.419580E-01
3   6.981317E-01   6.426908E-01   7.660918E-01 -6.427230E-01
4   1.396263E+00   9.847531E-01   1.737333E-01 -9.847453E-01
5   2.094395E+00 8.660529E-01 -4.998747E-01 -8.660057E-01
6   2.792527E+00   3.421830E-01 -9.395474E-01 -3.420648E-01
7   3.490659E+00 -3.417234E-01 -9.396111E-01   3.418948E-01
8   4.188790E+00 -8.656880E-01 -5.000588E-01   8.658733E-01
9   4.886922E+00 -9.845794E-01   1.734571E-01   9.847518E-01
10   5.585054E+00 -6.427721E-01   7.658258E-01   6.429526E-01
11   5.934120E+00 -3.420819E-01   9.395434E-01   3.423986E-01
12   6.283185E+00 -1.123186E-04   1.000000E+00   6.743190E-04
Error estimates     2.840430E-04   1.792939E-04   5.588399E-04

Additional Examples

Example 2

In this example, the following nonlinear problem is solved:

y″ − y3 + (1 + sin2t) sin t = 0

with y(0) = y(π) = 0. Its solution is y = sin t. As in Example 1, this equation is reduced to a system of first-order differential equations by defining y1 = y and y2= yʹ. The resulting system is

In this problem, there is one boundary condition at the left endpoint and one at the right endpoint; there are no coupled boundary conditions.

Note that since the parameter p is not used, in the call to BVPFD the routines FCNPEQ and FCNPBC are not needed. Therefore, in the call to BVPFD, FCNEQN and FCNBC were used in place of FCNPEQ and FCNPBC.

USE BVPFD_INT

USE UMACH_INT

USE CONST_INT

IMPLICIT NONE

! SPECIFICATIONS FOR PARAMETERS

INTEGER LDYFIN, LDYINI, MXGRID, NEQNS, NINIT

PARAMETER (MXGRID=45, NEQNS=2, NINIT=12, LDYFIN=NEQNS, &

LDYINI=NEQNS)

! SPECIFICATIONS FOR LOCAL VARIABLES

INTEGER I, J, NCUPBC, NFINAL, NLEFT, NOUT

REAL ERREST(NEQNS), PISTEP, TFINAL(MXGRID), TINIT(NINIT), &

TLEFT, TOL, TRIGHT, YFINAL(LDYFIN,MXGRID), &

YINIT(LDYINI,NINIT)

LOGICAL LINEAR, PRINT

! SPECIFICATIONS FOR INTRINSICS

INTRINSIC FLOAT

REAL FLOAT

! SPECIFICATIONS FOR FUNCTIONS

EXTERNAL FCNBC, FCNEQN, FCNJAC

! Set parameters

NLEFT = 1

NCUPBC = 0

TOL = .001

TLEFT = 0.0

TRIGHT = CONST('PI')

PISTEP = 0.0

PRINT = .FALSE.

LINEAR = .FALSE.

! Define TINIT and YINIT

DO 10 I=1, NINIT

TINIT(I) = TLEFT + (I-1)*(TRIGHT-TLEFT)/FLOAT(NINIT-1)

YINIT(1,I) = 0.4*(TINIT(I)-TLEFT)*(TRIGHT-TINIT(I))

YINIT(2,I) = 0.4*(TLEFT-TINIT(I)+TRIGHT-TINIT(I))

10 CONTINUE

! Solve problem

CALL BVPFD (FCNEQN, FCNJAC, FCNBC, FCNEQN, FCNBC, NLEFT, &

NCUPBC, TLEFT, TRIGHT, PISTEP, TOL, TINIT, &

YINIT, LINEAR, MXGRID, NFINAL, &

TFINAL, YFINAL, ERREST)

! Print results

CALL UMACH (2, NOUT)

WRITE (NOUT,99997)

WRITE (NOUT,99998) (I,TFINAL(I),(YFINAL(J,I),J=1,NEQNS),I=1, &

NFINAL)

WRITE (NOUT,99999) (ERREST(J),J=1,NEQNS)

99997 FORMAT (4X, 'I', 7X, 'T', 14X, 'Y1', 13X, 'Y2')

99998 FORMAT (I5, 1P3E15.6)

99999 FORMAT (' Error estimates', 4X, 1P2E15.6)

END

SUBROUTINE FCNEQN (NEQNS, T, Y, P, DYDT)

! SPECIFICATIONS FOR ARGUMENTS

INTEGER NEQNS

REAL T, P, Y(NEQNS), DYDT(NEQNS)

! SPECIFICATIONS FOR INTRINSICS

INTRINSIC SIN

REAL SIN

! Define PDE

DYDT(1) = Y(2)

DYDT(2) = Y(1)**3 - SIN(T)*(1.0+SIN(T)**2)

RETURN

END

SUBROUTINE FCNJAC (NEQNS, T, Y, P, DYPDY)

! SPECIFICATIONS FOR ARGUMENTS

INTEGER NEQNS

REAL T, P, Y(NEQNS), DYPDY(NEQNS,NEQNS)

! Define d(DYDT)/dY

DYPDY(1,1) = 0.0

DYPDY(1,2) = 1.0

DYPDY(2,1) = 3.0*Y(1)**2

DYPDY(2,2) = 0.0

RETURN

END

SUBROUTINE FCNBC (NEQNS, YLEFT, YRIGHT, P, F)

! SPECIFICATIONS FOR ARGUMENTS

INTEGER NEQNS

REAL P, YLEFT(NEQNS), YRIGHT(NEQNS), F(NEQNS)

! Define boundary conditions

F(1) = YLEFT(1)

F(2) = YRIGHT(1)

RETURN

END

Output

I       T              Y1             Y2
1   0.000000E+00   0.000000E+00   9.999277E-01
2   2.855994E-01   2.817682E-01   9.594315E-01
3   5.711987E-01   5.406458E-01   8.412407E-01
4   8.567980E-01   7.557380E-01   6.548904E-01
5   1.142397E+00   9.096186E-01   4.154530E-01
6   1.427997E+00   9.898143E-01   1.423307E-01
7   1.713596E+00   9.898143E-01 -1.423307E-01
8   1.999195E+00   9.096185E-01 -4.154530E-01
9   2.284795E+00   7.557380E-01 -6.548903E-01
10   2.570394E+00   5.406460E-01 -8.412405E-01
11   2.855994E+00   2.817683E-01 -9.594313E-01
12   3.141593E+00   0.000000E+00 -9.999274E-01
Error estimates     3.906105E-05   7.124186E-05

Example 3

In this example, the following nonlinear problem is solved:

with y(0) = y(1) = π/2. As in the previous examples, this equation is reduced to a system of first-order differential equations by defining y1 = y and y2 = yʹ. The resulting system is

The problem is embedded in a family of problems by introducing the parameter p and by changing the second differential equation to

At p = 0, the problem is linear; and at p = 1, the original problem is recovered. The derivatives
∂yʹ/∂p must now be specified in the subroutine FCNPEQ. The derivatives ∂f/∂p are zero in FCNPBC.

USE BVPFD_INT

USE UMACH_INT

IMPLICIT NONE

! SPECIFICATIONS FOR PARAMETERS

INTEGER LDYFIN, LDYINI, MXGRID, NEQNS, NINIT

PARAMETER (MXGRID=45, NEQNS=2, NINIT=5, LDYFIN=NEQNS, &

LDYINI=NEQNS)

! SPECIFICATIONS FOR LOCAL VARIABLES

INTEGER NCUPBC, NFINAL, NLEFT, NOUT

REAL ERREST(NEQNS), PISTEP, TFINAL(MXGRID), TLEFT, TOL, &

XRIGHT, YFINAL(LDYFIN,MXGRID)

LOGICAL LINEAR, PRINT

! SPECIFICATIONS FOR SAVE VARIABLES

INTEGER I, J

REAL TINIT(NINIT), YINIT(LDYINI,NINIT)

SAVE I, J, TINIT, YINIT

! SPECIFICATIONS FOR FUNCTIONS

EXTERNAL FCNBC, FCNEQN, FCNJAC, FCNPBC, FCNPEQ

DATA TINIT/0.0, 0.4, 0.5, 0.6, 1.0/

DATA ((YINIT(I,J),J=1,NINIT),I=1,NEQNS)/0.15749, 0.00215, 0.0, &

0.00215, 0.15749, -0.83995, -0.05745, 0.0, 0.05745, 0.83995/

! Set parameters

NLEFT = 1

NCUPBC = 0

TOL = .001

TLEFT = 0.0

XRIGHT = 1.0

PISTEP = 0.1

PRINT = .FALSE.

LINEAR = .FALSE.

CALL BVPFD (FCNEQN, FCNJAC, FCNBC, FCNPEQ, FCNPBC, NLEFT, &

NCUPBC, TLEFT, XRIGHT, PISTEP, TOL, TINIT, &

YINIT, LINEAR, MXGRID, NFINAL,TFINAL, YFINAL, ERREST)

! Print results

CALL UMACH (2, NOUT)

WRITE (NOUT,99997)

WRITE (NOUT,99998) (I,TFINAL(I),(YFINAL(J,I),J=1,NEQNS),I=1, &

NFINAL)

WRITE (NOUT,99999) (ERREST(J),J=1,NEQNS)

99997 FORMAT (4X, 'I', 7X, 'T', 14X, 'Y1', 13X, 'Y2')

99998 FORMAT (I5, 1P3E15.6)

99999 FORMAT (' Error estimates', 4X, 1P2E15.6)

END

SUBROUTINE FCNEQN (NEQNS, T, Y, P, DYDT)

! SPECIFICATIONS FOR ARGUMENTS

INTEGER NEQNS

REAL T, P, Y(NEQNS), DYDT(NEQNS)

! Define PDE

DYDT(1) = Y(2)

DYDT(2) = P*Y(1)**3 + 40./9.*((T-0.5)**2)**(1./3.) - (T-0.5)**8

RETURN

END

SUBROUTINE FCNJAC (NEQNS, T, Y, P, DYPDY)

! SPECIFICATIONS FOR ARGUMENTS

INTEGER NEQNS

REAL T, P, Y(NEQNS), DYPDY(NEQNS,NEQNS)

! Define d(DYDT)/dY

DYPDY(1,1) = 0.0

DYPDY(1,2) = 1.0

DYPDY(2,1) = P*3.*Y(1)**2

DYPDY(2,2) = 0.0

RETURN

END

SUBROUTINE FCNBC (NEQNS, YLEFT, YRIGHT, P, F)

USE CONST_INT

! SPECIFICATIONS FOR ARGUMENTS

INTEGER NEQNS

REAL P, YLEFT(NEQNS), YRIGHT(NEQNS), F(NEQNS)

! SPECIFICATIONS FOR LOCAL VARIABLES

REAL PI

! Define boundary conditions

PI = CONST('PI')

F(1) = YLEFT(1) - PI/2.0

F(2) = YRIGHT(1) - PI/2.0

RETURN

END

SUBROUTINE FCNPEQ (NEQNS, T, Y, P, DYPDP)

! SPECIFICATIONS FOR ARGUMENTS

INTEGER NEQNS

REAL T, P, Y(NEQNS), DYPDP(NEQNS)

! Define d(DYDT)/dP

DYPDP(1) = 0.0

DYPDP(2) = Y(1)**3

RETURN

END

SUBROUTINE FCNPBC (NEQNS, YLEFT, YRIGHT, P, DFDP)

! SPECIFICATIONS FOR ARGUMENTS

INTEGER NEQNS

REAL P, YLEFT(NEQNS), YRIGHT(NEQNS), DFDP(NEQNS)

! SPECIFICATIONS FOR SUBROUTINES

EXTERNAL SSET

! Define dF/dP

CALL SSET (NEQNS, 0.0, DFDP, 1)

RETURN

END

Output

I       T              Y1             Y2
1   0.000000E+00   1.570796E+00 -1.949336E+00
2   4.444445E-02   1.490495E+00 -1.669567E+00
3   8.888889E-02   1.421951E+00 -1.419465E+00
4   1.333333E-01   1.363953E+00 -1.194307E+00
5   2.000000E-01   1.294526E+00 -8.958461E-01
6   2.666667E-01   1.243628E+00 -6.373191E-01
7   3.333334E-01   1.208785E+00 -4.135206E-01
8   4.000000E-01   1.187783E+00 -2.219351E-01
9   4.250000E-01   1.183038E+00 -1.584200E-01
10   4.500000E-01   1.179822E+00 -9.973146E-02
11   4.625000E-01   1.178748E+00 -7.233893E-02
12   4.750000E-01   1.178007E+00 -4.638248E-02
13   4.812500E-01   1.177756E+00 -3.399763E-02
14   4.875000E-01   1.177582E+00 -2.205547E-02
15   4.937500E-01   1.177480E+00 -1.061177E-02
16   5.000000E-01   1.177447E+00 -1.479182E-07
17   5.062500E-01   1.177480E+00   1.061153E-02
18   5.125000E-01   1.177582E+00   2.205518E-02
19   5.187500E-01   1.177756E+00   3.399727E-02
20   5.250000E-01   1.178007E+00   4.638219E-02
21   5.375000E-01   1.178748E+00   7.233876E-02
22   5.500000E-01   1.179822E+00   9.973124E-02
23   5.750000E-01   1.183038E+00   1.584199E-01
24   6.000000E-01   1.187783E+00   2.219350E-01
25   6.666667E-01   1.208786E+00   4.135205E-01
26   7.333333E-01   1.243628E+00   6.373190E-01
27   8.000000E-01   1.294526E+00   8.958461E-01
28   8.666667E-01   1.363953E+00   1.194307E+00
29   9.111111E-01   1.421951E+00   1.419465E+00
30   9.555556E-01   1.490495E+00   1.669566E+00
31   1.000000E+00   1.570796E+00   1.949336E+00
Error estimates     3.448358E-06   5.549869E-05

Visual Numerics, Inc.
Visual Numerics - Developers of IMSL and PV-WAVE
http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260