## ODE: A Differential Equation Solver

``` /****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: ODE                                                 */
/*   TITLE: ODE: A Differential Equation Solver                 */
/* PRODUCT: IML                                                 */
/*  SYSTEM: ALL                                                 */
/*    KEYS:                                                     */
/*   PROCS: IML                                                 */
/* SUPPORT: GHG                         UPDATE:                 */
/*     REF:                                                     */
/*    MISC:                                                     */
/*                                                              */
/****************************************************************/

PROC IML  ;

/*-------------------------------------------------------------*/
/* EXAMPLE1: A SIMPLE PROGRAM TO SOLVE THE DIFFERENTIAL        */
/*           EQUATION                                          */
/*                                                             */
/*           dx                                                */
/*          --- = - t x     with x(0) = .5                     */
/*           dt                                                */
/*                                                             */
/* THE USER MUST SUPPLY                                        */
/*                                                             */
/*          1) THE (ORDERED) LIST OF POINTS AT WHICH THE       */
/*             SOLUTION IS TO BE  REPORTED.                    */
/*             THE INTEGRATION IS TO BE CARRIED UP TO .3       */
/*             AND ALSO REPORT THE RESULTS AT .1 and .2        */
/*                                                             */
/*          2) THE RIGHT HAND SIDE   (here it is : -t*x)       */
/*                                                             */
/*          3) THE STEP-SIZE VECTOR                            */
/*             THE INTEGRATION IS TO BE CARRIED  WITH          */
/*                 A MAXIMUM STEP SIZE OF 1                    */
/*                 A MINIMUM STEP SIZE OF 1E-12                */
/*                 AN INITIAL STEP SIZE OF 1E-5                */
/*          THESE NUMBERS ARE SPECIFIED IN H (NOT ORDERED)     */
/*                                                             */
/* THE USER CAN SUPPLY                                         */
/*                                                             */
/*          1) THE JACOBIAN OF THE RIGHT HAND SIDE  (HERE : -T)*/
/*          IF NOT SUPPLIED FINITE DIFFERENCE METHODS WILL     */
/*          BE  USED TO ESTIMATE THE JACOBIAN                  */
/*                                                             */
/* RESULTS ARE                                                 */
/*                                                             */
/*          Z[1,1] =  4.975062395981900000E-01                 */
/*          Z[1,2] =  4.900993366551400000E-01                 */
/*          Z[1,3] =  4.779987409188300000E-01                 */
/*-------------------------------------------------------------*/

START FUN(T,X) GLOBAL(COUNT);       /* DEFINE THE RHS          */
COUNT = COUNT+1;                /* WE COUNT THE NUMBER OF  */
/* FUNCTION EVALUATIONS    */
V = -T*X;
RETURN(V);
FINISH;

START JAC(T,X) ;                    /* DEFINE THE JACOBIAN     */
V = -T;
RETURN(V);
FINISH;

T = { 0 .1 .2 .3};                  /* SPECIFY THE ORDERED LIST*/
/* OF POINTS THAT DESCRIBE */
/* THE INTERVALS           */
C = .5;                             /* SPECIFY THE INITIAL     */
/* VALUE                   */

H = {1.E-12 1 1.E-5};               /* SPECIFY THE H VECTOR    */

COUNT = 0;                          /* SET THE COUNTER         */

CALL ODE(Z,"FUN",C,T,H) EPS = 1.E-12 J="JAC";
/* CALL ODE                */

PRINT COUNT;
PRINT Z[FORMAT =E25.15];

/*-------------------------------------------------------------*/
/* EXAMPLE2: A SIMPLE PROGRAM TO SOLVE THE DIFFERENTIAL        */
/*           EQUATION                                          */
/*                                                             */
/*           dx                        | 1 |                   */
/*          --- =  A x     with x(0) = |   |                   */
/*           dt                        | 2 |                   */
/*                                                             */
/*            |           |                                    */
/*            | -1      1 |                                    */
/*        A = |           |                                    */
/*            |  1      2 |                                    */
/*            |           |                                    */
/*                                                             */
/* THE USER MUST SUPPLY                                        */
/*                                                             */
/*          1) THE (ORDERED) LIST OF POINTS AT WHICH THE       */
/*             SOLUTION IS TO BE  REPORTED.                    */
/*             THE INTEGRATION IS TO BE CARRIED UP TO .3       */
/*             AND ALSO REPORT THE RESULTS AT .1 AND .2        */
/*                                                             */
/*          2) THE RIGHT HAND SIDE   (HERE IT IS : -T*X)       */
/*                                                             */
/*          3) THE STEP-SIZE VECTOR                            */
/*             THE INTEGRATION IS TO BE CARRIED  WITH          */
/*                 A MAXIMUM STEP SIZE OF 1                    */
/*                 A MINIMUM STEP SIZE OF 1E-12                */
/*                 AN INITIAL STEP SIZE OF 1E-5                */
/*          THESE NUMBERS ARE SPECIFIED IN H (NOT ORDERED)     */
/*                                                             */
/*          4) THE INITIAL VECTOR                              */
/*                                                             */
/* THE USER CAN SUPPLY                                         */
/*                                                             */
/*          1) THE JACOBIAN OF THE RIGHT HAND SIDE  (HERE : -T)*/
/*          IF NOT SUPPLIED FINITE DIFFERENCE METHODS WILL     */
/*          BE  USED TO ESTIMATE THE JACOBIAN                  */
/* RESULTS                                                     */
/*                                                             */
/*          FOR T= 1 THE RESULT IS STORED IN Z[,1]             */
/*                                                             */
/*            Z[1,1] = 1.059276279257100000E+00                */
/*            Z[2,1] = 7.554305542309500000E-01                */
/*                                                             */
/*          FOR T= 2 THE RESULT IS STORED IN Z[,2]             */
/*                                                             */
/*            Z[1,2] = 7.504499669077400000E-01                */
/*            Z[2,2] = 4.711431015044000000E-01                */
/*                                                             */
/*          FOR T= 3 THE RESULT IS STORED IN Z[,3]             */
/*                                                             */
/*            Z[1,3] = 5.141930081253500000E-01                */
/*            Z[2,3] = 3.183434053151600000E-01                */
/*                                                             */
/*-------------------------------------------------------------*/

START FUN2(T,X) GLOBAL(A) ;           /* DEFINE THE RHS        */
/* PASS A AS A PARAMETER */
V = A*X;
RETURN(V);
FINISH;

T = { 0  1  2  3 };                   /* DECLARE THE INTERVALS */

C = { 1 , 2 };                        /* DECLARE THE INTIAL    */
/* VECTOR C              */
A = { -1  1 , 1  -2 } ;               /* DEFINE A              */

H = {1.E-9 1 1.E-5};                  /* SET THE H VECTOR      */

CALL ODE(Z,"FUN2",C,T,H) EPS = 1.E-5; /* CALL THE INTEGRATOR   */

DO I = 1 TO 3;
D = Z[,I];
PRINT D[FORMAT =E25.15];
END;

```