fix import/export of ierode common on Windows
[scilab.git] / scilab / modules / differential_equations / src / fortran / rkf45.f
1       subroutine rkf45(fydot,neqn,y,t,tout,itol,rerr,aerr,
2      1    itask,iflag,iopt,work,lrw,iwork,liw,bjac,mf)
3 c
4 c     fehlberg fourth-fifth order runge-kutta method
5 c
6 c     written by h.a.watts and l.f.shampine
7 c                   sandia laboratories
8 c                  albuquerque,new mexico
9 c
10 c    rkf45 is primarily designed to solve non-stiff and mildly stiff
11 c    differential equations when derivative evaluations are inexpensive.
12 c    rkf45 should generally not be used when the user is demanding
13 c    high accuracy.
14 c
15 c abstract
16 c
17 c    subroutine  rkf45  integrates a system of neqn first order
18 c    ordinary differential equations of the form
19 c             dy(i)/dt = fydot(t,y(1),y(2),...,y(neqn))
20 c              where the y(i) are given at t .
21 c    typically the subroutine is used to integrate from t to tout but it
22 c    can be used as a one-step integrator to advance the solution a
23 c    single step in the direction of tout.  on return the parameters in
24 c    the call list are set for continuing the integration. the user has
25 c    only to call rkf45 again (and perhaps define a new value for tout).
26 c    actually, rkf45 is an interfacing routine which calls subroutine
27 c    rkfs for the solution.  rkfs in turn calls subroutine  fehl which
28 c    computes an approximate solution over one step.
29 c
30 c    rkf45  uses the runge-kutta-fehlberg (4,5)  method described
31 c    in the reference
32 c    e.fehlberg , low-order classical runge-kutta formulas with stepsize
33 c                 control , nasa tr r-315
34 c
35 c    the performance of rkf45 is illustrated in the reference
36 c    l.f.shampine,h.a.watts,s.davenport, solving non-stiff ordinary
37 c                 differential equations-the state of the art ,
38 c                 sandia laboratories report sand75-0182 ,
39 c                 to appear in siam review.
40 c
41 c
42 c    the parameters represent-
43 c      fydot -- subroutine fydot(neqn,t,y,yp) to evaluate derivatives 
44 c             yp(i)=dy(i)/dt
45 c      neqn -- number of equations to be integrated
46 c      y(*) -- solution vector at t
47 c      t -- independent variable
48 c      tout -- output point at which solution is desired
49 c      rerr,aerr -- relative and absolute error tolerances for local
50 c            error test. at each step the code requires that
51 c                 abs(local error) .le. rerr*abs(y) + aerr
52 c            for each component of the local error and solution vectors
53 c      iflag -- indicator for status of integration
54 c      work(*) -- array to hold information internal to rkf45 which is
55 c            necessary for subsequent calls. must be dimensioned
56 c            at least  3+6*neqn
57 c      iwork(*) -- integer array used to hold information internal to
58 c            rkf45 which is necessary for subsequent calls. must be
59 c            dimensioned at least  5
60 c
61 c
62 c  first call to rkf45
63 c
64 c    the user must provide storage in his calling program for the arrays
65 c    in the call list  -      y(neqn) , work(3+6*neqn) , iwork(5)  ,
66 c    declare fydot in an external statement, supply 
67 c    subroutine fydot(neqn,t,y,yp)
68 c    and initialize the following parameters-
69 c
70 c      neqn -- number of equations to be integrated.  (neqn .ge. 1)
71 c      y(*) -- vector of initial conditions
72 c      t -- starting point of integration , must be a variable
73 c      tout -- output point at which solution is desired.
74 c            t=tout is allowed on the first call only, in which case
75 c            rkf45 returns with iflag=2 if continuation is possible.
76 c      rerr,aerr -- relative and absolute local error tolerances
77 c            which must be non-negative. rerr must be a variable while
78 c            aerr may be a constant. the code should normally not be
79 c            used with relative error control smaller than about 1.e-8 .
80 c            to avoid limiting precision difficulties the code requires
81 c            rerr to be larger than an internally computed relative
82 c            error parameter which is machine dependent. in particular,
83 c            pure absolute error is not permitted. if a smaller than
84 c            allowable value of rerr is attempted, rkf45 increases
85 c            rerr appropriately and returns control to the user before
86 c            continuing the integration.
87 c      iflag -- +1,-1  indicator to initialize the code for each new
88 c            problem. normal input is +1. the user should set iflag=-1
89 c            only when one-step integrator control is essential. in this
90 c            case, rkf45 attempts to advance the solution a single step
91 c            in the direction of tout each time it is called. since this
92 c            mode of operation results in extra computing overhead, it
93 c            should be avoided unless needed.
94 c
95 c
96 c  output from rkf45
97 c
98 c      y(*) -- solution at t
99 c      t -- last point reached in integration.
100 c      iflag = 2 -- integration reached tout. indicates successful retur
101 c                   and is the normal mode for continuing integration.
102 c            =-2 -- a single successful step in the direction of tout
103 c                   has been taken. normal mode for continuing
104 c                   integration one step at a time.
105 c            = 3 -- integration was not completed because relative error
106 c                   tolerance was too small. rerr has been increased
107 c                   appropriately for continuing.
108 c            = 4 -- integration was not completed because more than
109 c                   3000 derivative evaluations were needed. this
110 c                   is approximately 500 steps.
111 c            = 5 -- integration was not completed because solution
112 c                   vanished making a pure relative error test
113 c                   impossible. must use non-zero aerr to continue.
114 c                   using the one-step integration mode for one step
115 c                   is a good way to proceed.
116 c            = 6 -- integration was not completed because requested
117 c                   accuracy could not be achieved using smallest
118 c                   allowable stepsize. user must increase the error
119 c                   tolerance before continued integration can be
120 c                   attempted.
121 c            = 7 -- it is likely that rkf45 is inefficient for solving
122 c                   this problem. too much output is restricting the
123 c                   natural stepsize choice. use the one-step integrator
124 c                   mode.
125 c            = 8 -- invalid input parameters
126 c                   this indicator occurs if any of the following is
127 c                   satisfied -   neqn .le. 0
128 c                                 t=tout  and  iflag .ne. +1 or -1
129 c                                 rerr or aerr .lt. 0.
130 c                                 iflag .eq. 0  or  .lt. -2  or  .gt. 8
131 c      work(*),iwork(*) -- information which is usually of no interest
132 c                   to the user but necessary for subsequent calls.
133 c                   work(1),...,work(neqn) contain the first derivatives
134 c                   of the solution vector y at t. work(neqn+1) contains
135 c                   the stepsize h to be attempted on the next step.
136 c                   iwork(1) contains the derivative evaluation counter.
137 c
138 c
139 c  subsequent calls to rkf45
140 c
141 c    subroutine rkf45 returns with all information needed to continue
142 c    the integration. if the integration reached tout, the user only
143 c    needs to define a new tout and call rkf45 again. in the one-step
144 c    integrator mode (iflag=-2) the user must keep in mind that each
145 c    step taken is in the direction of the current tout. upon reaching
146 c    tout (indicated by changing iflag to 2),the user must then define
147 c    a new tout and reset iflag to -2 to continue in the one-step
148 c    integrator mode.
149 c
150 c    if the integration was not completed but the user still wants to
151 c    continue (iflag=3,4 cases), he just calls rkf45 again. with iflag=3
152 c    the rerr parameter has been adjusted appropriately for continuing
153 c    the integration. in the case of iflag=4 the function counter will
154 c    be reset to 0 and another 3000 function evaluations are allowed.
155 c
156 c    however,in the case iflag=5, the user must first alter the error
157 c    criterion to use a positive value of aerr before integration can
158 c    proceed. if he does not,execution is terminated.
159 c
160 c    also,in the case iflag=6, it is necessary for the user to reset
161 c    iflag to 2 (or -2 when the one-step integration mode is being used)
162 c    as well as increasing either aerr,rerr or both before the
163 c    integration can be continued. if this is not done, execution will
164 c    be terminated. the occurrence of iflag=6 indicates a trouble spot
165 c    (solution is changing rapidly,singularity may be present) and it
166 c    often is inadvisable to continue.
167 c
168 c    if iflag=7 is encountered, the user should use the one-step
169 c    integration mode with the stepsize determined by the code or
170 c    consider switching to the adams codes de/step,intrp. if the user
171 c    insists upon continuing the integration with rkf45, he must reset
172 c    iflag to 2 before calling rkf45 again. otherwise,execution will be
173 c    terminated.
174 c
175 c    if iflag=8 is obtained, integration can not be continued unless
176 c    the invalid input parameters are corrected.
177 c
178 c    it should be noted that the arrays work,iwork contain information
179 c    required for subsequent integration. accordingly, work and iwork
180 c    should not be altered.
181 c
182 c
183       include 'stack.h'
184       
185       integer neqn,iflag,iwork(5)
186       double precision y(neqn),t,tout,rerr,aerr,work(1)
187 c
188       external fydot
189 c
190       integer k1,k2,k3,k4,k5,k6,k1m
191       if(itask.eq.1) then iflag=1
192       if(itask.eq.2) then iflag=-1
193 c
194 c
195 c     compute indices for the splitting of the work array
196 c
197       k1m=neqn+1
198       k1=k1m+1
199       k2=k1+neqn
200       k3=k2+neqn
201       k4=k3+neqn
202       k5=k4+neqn
203       k6=k5+neqn
204       k7=k6+neqn
205 c
206 c     this interfacing routine merely relieves the user of a long
207 c     calling list via the splitting apart of two working storage
208 c     arrays. if this is not compatible with the users compiler,
209 c     he must use rkfs directly.
210 c
211       call rkfs(fydot,neqn,y,t,tout,rerr,aerr,iflag,work(1),work(k1m),
212      1          work(k1),work(k2),work(k3),work(k4),work(k5),work(k6),
213      2          work(k6+1),work(k7),
214      3          iwork(1),iwork(2),iwork(3),iwork(4),iwork(5))
215 c
216       return
217       end
218       subroutine rkfs(fydot,neqn,y,t,tout,rerr,aerr,iflag,yp,h,f1,f2,f3,
219      1                f4,f5,savre,savae,savey,nfe,kop,init,jflag,kflag)
220 c
221 c     fehlberg fourth-fifth order runge-kutta method
222 c
223 c
224 c     rkfs integrates a system of first order ordinary differential
225 c     equations as described in the comments for rkf45 .
226 c     the arrays yp,f1,f2,f3,f4,and f5 (of dimension at least neqn) and
227 c     the variables h,savre,savae,nfe,kop,init,jflag,and kflag are used
228 c     internally by the code and appear in the call list to eliminate
229 c     local retention of variables between calls. accordingly, they
230 c     should not be altered. items of possible interest are
231 c         yp - derivative of solution vector at t
232 c         h  - an appropriate stepsize to be used for the next step
233 c         nfe- counter on the number of derivative function evaluations
234 c
235 c
236       
237       logical hfaild,output
238 c
239       integer  neqn,iflag,nfe,kop,init,jflag,kflag
240       double precision  y(neqn),t,tout,rerr,aerr,h,yp(neqn),
241      1  f1(neqn),f2(neqn),f3(neqn),f4(neqn),f5(neqn),savre,
242      2  savae,savey(*)
243 c
244       external fydot
245 c
246       double precision  a,ae,dt,ee,eeoet,esttol,et,hmin,remin,rer,s,
247      1  scale,tol,toln,twoeps,u26,ypk
248 c
249       integer  k,maxnfe,mflag
250 c
251       double precision  dabs,dmax1,dmin1,dsign,dlamch
252 c
253 c  remin is the minimum acceptable value of rerr.  attempts
254 c  to obtain higher accuracy with this subroutine are usually
255 c  very expensive and often unsuccessful.
256 c
257       data remin/1.d-12/
258 c
259 c
260 c     the expense is controlled by restricting the number
261 c     of function evaluations to be approximately maxnfe.
262 c     as set, this corresponds to about 500 steps.
263 c
264       data maxnfe/3000/
265 c
266 c   here two constants emboding the machine epsilon is present
267 c   twoesp is set to twice the machine epsilon while u26 is set
268 c   to 26 times the machine epsilon
269 c
270 c     data twoeps, u26/4.4d-16, 5.72d-15/
271       twoeps = 2.*dlamch('p')
272       u26 = 13.*twoeps
273 c
274 c
275 c     check input parameters
276 c
277 c
278       if (neqn .lt. 1) go to 10
279       if ((rerr .lt. 0.0d0)  .or.  (aerr .lt. 0.0d0)) go to 10
280       mflag=iabs(iflag)
281       if ((mflag .ge. 1)  .and.  (mflag .le. 8)) go to 20
282 c
283 c     invalid input
284    10 iflag=8
285       return
286 c
287 c     is this the first call
288    20 if (mflag .eq. 1) go to 50
289 c
290 c     check continuation possibilities
291 c
292       if ((t .eq. tout) .and. (kflag .ne. 3)) go to 10
293       if (mflag .ne. 2) go to 25
294 c
295 c     iflag = +2 or -2
296       if (kflag .eq. 3) go to 45
297       if (init .eq. 0) go to 45
298       if (kflag .eq. 4) go to 40
299       if ((kflag .eq. 5)  .and.  (aerr .eq. 0.0d0)) go to 30
300       if ((kflag .eq. 6)  .and.  (rerr .le. savre)  .and.
301      1    (aerr .le. savae)) go to 30
302       go to 50
303 c
304 c     iflag = 3,4,5,6,7 or 8
305    25 if (iflag .eq. 3) go to 45
306       if (iflag .eq. 4) go to 40
307       if ((iflag .eq. 5) .and. (aerr .gt. 0.0d0)) go to 45
308 c
309 c     integration cannot be continued since user did not respond to
310 c     the instructions pertaining to iflag=5,6,7 or 8
311    30 stop
312 c
313 c     reset function evaluation counter
314    40 nfe=0
315       if (mflag .eq. 2) go to 50
316 c
317 c     reset flag value from previous call
318    45 iflag=jflag
319       if (kflag .eq. 3) mflag=iabs(iflag)
320 c
321 c     save input iflag and set continuation flag value for subsequent
322 c     input checking
323    50 jflag=iflag
324       kflag=0
325 c
326 c     save rerr and aerr for checking input on subsequent calls
327       savre=rerr
328       savae=aerr
329 c
330 c     restrict relative error tolerance to be at least as large as
331 c     2*eps+remin to avoid limiting precision difficulties arising
332 c     from impossible accuracy requests
333 c
334       rer=twoeps+remin
335       if (rerr .ge. rer) go to 55
336 c
337 c     relative error tolerance too small
338       rerr=rer
339       iflag=3
340       kflag=3
341       return
342 c
343    55 dt=tout-t
344 c
345       if (mflag .eq. 1) go to 60
346       if (init .eq. 0) go to 65
347       go to 80
348 c
349 c     initialization --
350 c                       set initialization completion indicator,init
351 c                       set indicator for too many output points,kop
352 c                       evaluate initial derivatives
353 c                       set counter for function evaluations,nfe
354 c                       evaluate initial derivatives
355 c                       set counter for function evaluations,nfe
356 c                       estimate starting stepsize
357 c
358    60 init=0
359       kop=0
360 c
361       a=t
362       call fydot(neqn,a,y,yp)
363       if(ierror.gt.0) return
364       nfe=1
365       if (t .ne. tout) go to 65
366       iflag=2
367       return
368 c
369 c
370    65 init=1
371       h=dabs(dt)
372       toln=0.
373       do 70 k=1,neqn
374         tol=rerr*dabs(y(k))+aerr
375         if (tol .le. 0.) go to 70
376         toln=tol
377         ypk=dabs(yp(k))
378         if (ypk*h**5 .gt. tol) h=(tol/ypk)**0.2d0
379    70 continue
380       if (toln .le. 0.0d0) h=0.0d0
381       h=dmax1(h,u26*dmax1(dabs(t),dabs(dt)))
382       jflag=isign(2,iflag)
383 c
384 c
385 c     set stepsize for integration in the direction from t to tout
386 c
387    80 h=dsign(h,dt)
388 c
389 c     test to see if rkf45 is being severely impacted by too many
390 c     output points
391 c
392       if (dabs(h) .ge. 2.0d0*dabs(dt)) kop=kop+1
393       if (kop .ne. 100) go to 85
394 c
395 c     unnecessary frequency of output
396       kop=0
397       iflag=7
398       return
399 c
400    85 if (dabs(dt) .gt. u26*dabs(t)) go to 95
401 c
402 c     if too close to output point,extrapolate and return
403 c
404       do 90 k=1,neqn
405    90   y(k)=y(k)+dt*yp(k)
406       a=tout
407       call fydot(neqn,a,y,yp)
408       if(ierror.gt.0) return
409       nfe=nfe+1
410       go to 300
411 c
412 c
413 c     initialize output point indicator
414 c
415    95 output= .false.
416 c
417 c     to avoid premature underflow in the error tolerance function,
418 c     scale the error tolerances
419 c
420       scale=2.0d0/rerr
421       ae=scale*aerr
422 c
423 c
424 c     step by step integration
425 c
426   100 hfaild= .false.
427 c
428 c     set smallest allowable stepsize
429 c
430       hmin=u26*dabs(t)
431 c
432 c     adjust stepsize if necessary to hit the output point.
433 c     look ahead two steps to avoid drastic changes in the stepsize and
434 c     thus lessen the impact of output points on the code.
435 c
436       dt=tout-t
437       if (dabs(dt) .ge. 2.0d0*dabs(h)) go to 200
438       if (dabs(dt) .gt. dabs(h)) go to 150
439 c
440 c     the next successful step will complete the integration to the
441 c     output point
442 c
443       output= .true.
444       h=dt
445       go to 200
446 c
447   150 h=0.5d0*dt
448 c
449 c
450 c
451 c     core integrator for taking a single step
452 c
453 c     the tolerances have been scaled to avoid premature underflow in
454 c     computing the error tolerance function et.
455 c     to avoid problems with zero crossings,relative error is measured
456 c     using the average of the magnitudes of the solution at the
457 c     beginning and end of a step.
458 c     the error estimate formula has been grouped to control loss of
459 c     significance.
460 c     to distinguish the various arguments, h is not permitted
461 c     to become smaller than 26 units of roundoff in t.
462 c     practical limits on the change in the stepsize are enforced to
463 c     smooth the stepsize selection process and to avoid excessive
464 c     chattering on problems having discontinuities.
465 c     to prevent unnecessary failures, the code uses 9/10 the stepsize
466 c     it estimates will succeed.
467 c     after a step failure, the stepsize is not allowed to increase for
468 c     the next attempted step. this makes the code more efficient on
469 c     problems having discontinuities and more effective in general
470 c     since local extrapolation is being used and extra caution seems
471 c     warranted.
472 c
473 c
474 c     test number of derivative function evaluations.
475 c     if okay,try to advance the integration from t to t+h
476 c
477   200 if (nfe .le. maxnfe) go to 220
478 c
479 c     too much work
480       iflag=4
481       kflag=4
482       return
483 c
484 c     advance an approximate solution over one step of length h
485 c
486   220 continue
487       do 33 k=1,neqn
488  33   savey(k)=y(k)
489       call fehl(fydot,neqn,y,t,h,yp,f1,f2,f3,f4,f5,f1,savey)
490       do 34 k=1,neqn
491  34   y(k)=savey(k)
492       nfe=nfe+5
493 c
494 c     compute and test allowable tolerances versus local error estimates
495 c     and remove scaling of tolerances. note that relative error is
496 c     measured with respect to the average of the magnitudes of the
497 c     solution at the beginning and end of the step.
498 c
499       eeoet=0.0d0
500       do 250 k=1,neqn
501         et=dabs(savey(k))+dabs(f1(k))+ae
502         if (et .gt. 0.0d0) go to 240
503 c
504 c       inappropriate error tolerance
505         iflag=5
506         return
507 c
508   240   ee=dabs((-2090.0d0*yp(k)+(21970.0d0*f3(k)-15048.0d0*f4(k)))+
509      1                        (22528.0d0*f2(k)-27360.0d0*f5(k)))
510   250   eeoet=dmax1(eeoet,ee/et)
511 c
512       esttol=dabs(h)*eeoet*scale/752400.0d0
513 c
514       if (esttol .le. 1.0d0) go to 260
515 c
516 c
517 c     unsuccessful step
518 c                       reduce the stepsize , try again
519 c                       the decrease is limited to a factor of 1/10
520 c
521       hfaild= .true.
522       output= .false.
523       s=0.1d0
524       if (esttol .lt. 59049.0d0) s=0.9d0/esttol**0.2d0
525       h=s*h
526       if (dabs(h) .gt. hmin) go to 200
527 c
528 c     requested error unattainable at smallest allowable stepsize
529       iflag=6
530       kflag=6
531       return
532 c
533 c
534 c     successful step
535 c                        store solution at t+h
536 c                        and evaluate derivatives there
537 c
538   260 t=t+h
539       do 270 k=1,neqn
540   270   y(k)=f1(k)
541       a=t
542       call fydot(neqn,a,y,yp)
543       if(ierror.gt.0) return
544       nfe=nfe+1
545 c
546 c
547 c                       choose next stepsize
548 c                       the increase is limited to a factor of 5
549 c                       if step failure has just occurred, next
550 c                          stepsize is not allowed to increase
551 c
552       s=5.0d0
553       if (esttol .gt. 1.889568d-4) s=0.9d0/esttol**0.2d0
554       if (hfaild) s=dmin1(s,1.0d0)
555       h=dsign(dmax1(s*dabs(h),hmin),h)
556 c
557 c     end of core integrator
558 c
559 c
560 c     should we take another step
561 c
562       if (output) go to 300
563       if (iflag .gt. 0) go to 100
564 c
565 c
566 c     integration successfully completed
567 c
568 c     one-step mode
569       iflag=-2
570       return
571 c
572 c     interval mode
573   300 t=tout
574       iflag=2
575       return
576 c
577       end
578       subroutine fehl(fydot,neqn,y,t,h,yp,f1,f2,f3,f4,f5,s,savey)
579 c      subroutine fehl(fydot,neqn,y,t,h,yp,f1,f2,f3,f4,f5,s)
580 c
581 c     fehlberg fourth-fifth order runge-kutta method
582 c
583 c    fehl integrates a system of neqn first order
584 c    ordinary differential equations of the form
585 c             dy(i)/dt=fydot(t,y(1),---,y(neqn))
586 c    where the initial values y(i) and the initial derivatives
587 c    yp(i) are specified at the starting point t. fehl advances
588 c    the solution over the fixed step h and returns
589 c    the fifth order (sixth order accurate locally) solution
590 c    approximation at t+h in array s(i).
591 c    f1,---,f5 are arrays of dimension neqn which are needed
592 c    for internal storage.
593 c    the formulas have been grouped to control loss of significance.
594 c    fehl should be called with an h not smaller than 13 units of
595 c    roundoff in t so that the various independent arguments can be
596 c    distinguished.
597 c
598 c
599       integer  neqn
600       double precision  y(neqn),t,h,yp(neqn),f1(neqn),f2(neqn),
601      1  f3(neqn),f4(neqn),f5(neqn),s(neqn),savey(neqn)
602 c
603       
604       include 'stack.h'
605       
606       double precision  ch
607       integer  k
608       external fydot
609 c
610       ch=h/4.0d0
611       do 221 k=1,neqn
612   221   y(k)=savey(k)+ch*yp(k)
613       call fydot(neqn,t+ch,y,f1)
614       if(ierror.gt.0) return
615 c
616       ch=3.0d0*h/32.0d0
617       do 222 k=1,neqn
618   222   y(k)=savey(k)+ch*(yp(k)+3.0d0*f1(k))
619       call fydot(neqn,t+3.0d0*h/8.0d0,y,f2)
620       if(ierror.gt.0) return
621 c
622       ch=h/2197.0d0
623       do 223 k=1,neqn
624   223 y(k)=savey(k)+ch*(1932.0d0*yp(k)+
625      1    (7296.0d0*f2(k)-7200.0d0*f1(k)))
626       call fydot(neqn,t+12.0d0*h/13.0d0,y,f3)
627       if(ierror.gt.0) return
628 c
629       ch=h/4104.0d0
630       do 224 k=1,neqn
631   224   y(k)=savey(k)+ch*((8341.0d0*yp(k)-845.0d0*f3(k))+
632      1                            (29440.0d0*f2(k)-32832.0d0*f1(k)))
633       call fydot(neqn,t+h,y,f4)
634       if(ierror.gt.0) return
635 c
636       ch=h/20520.0d0
637       do 225 k=1,neqn
638   225   y(k)=savey(k)+ch*((-6080.0d0*yp(k)+(9295.0d0*f3(k)-
639      1         5643.0d0*f4(k)))+(41040.0d0*f1(k)-28352.0d0*f2(k)))
640       call fydot(neqn,t+h/2.0d0,y,f5)
641       if(ierror.gt.0) return
642 c
643 c     compute approximate solution at t+h
644 c
645       ch=h/7618050.0d0
646       do 230 k=1,neqn
647   230   s(k)=savey(k)+ch*((902880.0d0*yp(k)+(3855735.0d0*f3(k)-
648      1        1371249.0d0*f4(k)))+(3953664.0d0*f2(k)+
649      2        277020.0d0*f5(k)))
650 c
651       return
652       end
653