#include <S.h>
extern void F77_NAME(heun)() ;
static void *Sdydt ; /* pointer to S-PLUS function to be filled in */
/* descriptions of the functions's two arguments */
static char *modes[] = { "single", "single" } ;
static long lengths[] = { 1, 0 } ; /* neqn = lengths[1] to be filled in */
static char *names[] = { "t", "y" } ;

static void
f(t, y, yp)
float *t ; /* 1 long, input */
float *y ; /* neqn long, input */
float *yp ; /* neqn long, output */
{
        void *in[2] ; /* for two inputs to S-PLUS function, t and y */
        void *out[1] ; /* for one output vector of S-PLUS function */
        int i;
        in[0] = (void *)t ;
        in[1] = (void *)y ;
        call_S(Sdydt,
                2L, /* 2 arguments */
                in, modes, lengths, names, /* all of length 2, 1 for each arg */
                1L, /* 1 result */
                out) ; /* the return value */
        /*
         * out must be 1 long---i.e., S-PLUS function must return
         * atomic vector or a list of one atomic vector.  We can
         * check that it is at least 1 long.
         */
        if (!out[0])
                PROBLEM
                        "S-PLUS function returned 0 long list"
                RECOVER(NULL_ENTRY) ;
        /*
         * Assume out[0] points to lengths[1] single precision numbers.
         * We cannot check this assumption here.
         */
        for(i=0;i<lengths[1];i++)
                yp[i] = ((float *)out[0])[i] ;
        return ;
}

/* called via .C by S-PLUS function dfeq() */
dfeq(Sdydt_p, y, neqn, t_start, t_end, step, work)
void **Sdydt_p ;
float *y;
long *neqn ;
float *t_start, *t_end, *step, *work ;
{
        /* Store pointer to S-PLUS function and number of equations */
        Sdydt = *Sdydt_p ;
        lengths[1] = *neqn ;
        /* call Fortran differential equation solver */
        F77_CALL(heun)(f, neqn, y, t_start, t_end, step, work) ;
}