b746f732e6301a48270e530df7148ce0e82ec879
[scilab.git] / scilab / modules / differential_equations / src / fortran / odeint.f
1 c     ====================================
2       subroutine odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,rk
3      *qc)
4 C      implicit undefined (a-z)
5       external derivs,rkqc
6       integer maxstp,nmax,kmax,kount,nvar,i,nok,nbad,nstp
7       double precision two,zero,tiny,dxsav,xp(200),yp(10,200),x,h
8       parameter (maxstp=10000,nmax=10,two=2.0,zero=0.0,tiny=1.e-30)
9       double precision x1,x2,eps,h1,hmin,xsav,hdid,hnext
10       double precision ystart(nvar),yscal(nmax),y(nmax),dydx(nmax)
11       character*80 messag
12       common /path/ kmax,kount,dxsav,xp,yp
13       integer iero
14       common/ierode/iero
15       iero=0
16       if ( abs(x2-x1).le.tiny) return
17       x=x1
18       h=sign(h1,x2-x1)
19       nok=0
20       nbad=0
21       kount=0
22       do 11 i=1,nvar
23         y(i)=ystart(i)
24 11    continue
25       xsav=x-dxsav*two
26       do 16 nstp=1,maxstp
27         call derivs(nvar,x,y,dydx)
28         if (iero.gt.0) return 
29         do 12 i=1,nvar
30           yscal(i)=abs(y(i))+abs(h*dydx(i))+tiny
31 12      continue
32         if(kmax.gt.0)then
33           if(abs(x-xsav).gt.abs(dxsav)) then
34             if(kount.lt.kmax-1)then
35               kount=kount+1
36               xp(kount)=x
37               do 13 i=1,nvar
38                 yp(i,kount)=y(i)
39 13            continue
40               xsav=x
41             endif
42           endif
43         endif
44         if((x+h-x2)*(x+h-x1).gt.zero) h=x2-x
45         call rkqc(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs)
46         if(iero.gt.0) return
47         if(hdid.eq.h)then
48           nok=nok+1
49         else
50           nbad=nbad+1
51         endif
52         if((x-x2)*(x2-x1).ge.zero)then
53           do 14 i=1,nvar
54             ystart(i)=y(i)
55 14        continue
56           if(kmax.ne.0)then
57             kount=kount+1
58             xp(kount)=x
59             do 15 i=1,nvar
60               yp(i,kount)=y(i)
61 15          continue
62           endif
63           return
64         endif
65         if(abs(hnext).lt.hmin) then
66            write(messag, 17) hnext,hmin
67            hnext=hmin
68         endif
69         h=hnext
70 16    continue
71       iero=-1
72 c      print *, 'Trop d''iterations a faire pour la precision demandee.'
73       return
74  17   format('stepsize ',e10.3,' smaller than minimum ',e10.3)
75       end
76 c     ====================================
77       
78