1 SUBROUTINE sstode (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR,
2 1 wm, iwm,
f, jac, pjac, slvs)
95 EXTERNAL f, jac, pjac, slvs
97 REAL Y, YH, YH1, EWT, SAVF, ACOR, WM
98 dimension neq(*), y(*), yh(nyh,*), yh1(*), ewt(*), savf(*),
99 1 acor(*), wm(*), iwm(*)
100 INTEGER INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, CNYH,
101 1 ialth, ipup, lmax, meo, nqnyh, nslp,
102 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
103 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
104 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
105 INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ
106 REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
107 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
108 REAL DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP,
109 1 r, rh, rhdn, rhsm, rhup, told, svnorm
110 COMMON /sls001/ conit, crate, el(13), elco(13,12),
111 1 hold, rmax, tesco(3,12),
112 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
113 2 init, mxstep, mxhnil, nhnil, nslast, cnyh,
114 3 ialth, ipup, lmax, meo, nqnyh, nslp,
115 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
116 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
117 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
128 IF (jstart .GT. 0) go
to 200
129 IF (jstart .EQ. -1) go
to 100
130 IF (jstart .EQ. -2) go
to 160
168 IF (ialth .EQ. 1) ialth = 2
169 IF (meth .EQ. meo) go
to 110
170 CALL
scfode(meth, elco, tesco)
172 IF (nq .GT. maxord) go
to 120
176 110
IF (nq .LE. maxord) go
to 160
180 125 el(i) = elco(i,nq)
185 ddn = svnorm(n, savf, ewt)/tesco(1,l)
187 rhdn = 1.0e0/(1.3e0*ddn**exdn + 0.0000013e0)
190 IF (h .EQ. hold) go
to 170
199 140 CALL
scfode(meth, elco, tesco)
201 155 el(i) = elco(i,nq)
206 go
to(160, 170, 200), iret
213 160
IF (h .EQ. hold) go
to 200
218 170 rh =
max(rh,hmin/
abs(h))
219 175 rh =
min(rh,rmax)
220 rh = rh/
max(1.0e0,
abs(h)*hmxi*rh)
225 180 yh(i,j) = yh(i,j)*r
229 IF (iredo .EQ. 0) go
to 690
238 200
IF (
abs(rc-1.0e0) .GT. ccmax) ipup = miter
239 IF (nst .GE. nslp+msbp) ipup = miter
246 210 yh1(i) = yh1(i) + yh1(i+nyh)
257 CALL
f(neq, tn, y, savf)
259 IF (ipup .LE. 0) go
to 250
265 CALL pjac(neq, y, yh, nyh, ewt, acor, savf, wm, iwm,
f, jac)
270 IF (ierpj .NE. 0) go
to 430
273 270
IF (miter .NE. 0) go
to 350
279 savf(i) = h*savf(i) - yh(i,2)
280 290 y(i) = savf(i) - acor(i)
281 del = svnorm(n, y, ewt)
283 y(i) = yh(i,1) + el(1)*savf(i)
284 300 acor(i) = savf(i)
292 360 y(i) = h*savf(i) - (yh(i,2) + acor(i))
293 CALL slvs(wm, iwm, y, savf)
294 IF (iersl .LT. 0) go
to 430
295 IF (iersl .GT. 0) go
to 410
296 del = svnorm(n, y, ewt)
298 acor(i) = acor(i) + y(i)
299 380 y(i) = yh(i,1) + el(1)*acor(i)
304 400
IF (m .NE. 0) crate =
max(0.2e0*crate,del/delp)
305 dcon = del*
min(1.0e0,1.5e0*crate)/(tesco(2,nq)*conit)
306 IF (dcon .LE. 1.0e0) go
to 450
308 IF (m .EQ. maxcor) go
to 410
309 IF (m .GE. 2 .AND. del .GT. 2.0e0*delp) go
to 410
311 CALL
f(neq, tn, y, savf)
321 410
IF (miter .EQ. 0 .OR. jcur .EQ. 1) go
to 430
334 440 yh1(i) = yh1(i) - yh1(i+nyh)
336 IF (ierpj .LT. 0 .OR. iersl .LT. 0) go
to 680
337 IF (
abs(h) .LE. hmin*1.00001e0) go
to 670
338 IF (ncf .EQ. mxncf) go
to 670
350 IF (m .EQ. 0) dsm = del/tesco(2,nq)
351 IF (m .GT. 0) dsm = svnorm(n, acor, ewt)/tesco(2,nq)
352 IF (dsm .GT. 1.0e0) go
to 500
370 470 yh(i,j) = yh(i,j) + el(j)*acor(i)
372 IF (ialth .EQ. 0) go
to 520
373 IF (ialth .GT. 1) go
to 700
374 IF (l .EQ. lmax) go
to 700
376 490 yh(i,lmax) = acor(i)
385 500 kflag = kflag - 1
392 510 yh1(i) = yh1(i) - yh1(i+nyh)
395 IF (
abs(h) .LE. hmin*1.00001e0) go
to 660
396 IF (kflag .LE. -3) go
to 640
410 IF (l .EQ. lmax) go
to 540
412 530 savf(i) = acor(i) - yh(i,lmax)
413 dup = svnorm(n, savf, ewt)/tesco(3,nq)
415 rhup = 1.0e0/(1.4e0*dup**exup + 0.0000014e0)
417 rhsm = 1.0e0/(1.2e0*dsm**exsm + 0.0000012e0)
419 IF (nq .EQ. 1) go
to 560
420 ddn = svnorm(n, yh(1,l), ewt)/tesco(1,nq)
422 rhdn = 1.0e0/(1.3e0*ddn**exdn + 0.0000013e0)
423 560
IF (rhsm .GE. rhup) go
to 570
424 IF (rhup .GT. rhdn) go
to 590
426 570
IF (rhsm .LT. rhdn) go
to 580
432 IF (kflag .LT. 0 .AND. rh .GT. 1.0e0) rh = 1.0e0
436 IF (rh .LT. 1.1e0) go
to 610
439 600 yh(i,newq+1) = acor(i)*r
443 620
IF ((kflag .EQ. 0) .AND. (rh .LT. 1.1e0)) go
to 610
444 IF (kflag .LE. -2) rh =
min(rh,0.2e0)
450 IF (newq .EQ. nq) go
to 170
464 640
IF (kflag .EQ. -10) go
to 660
470 CALL
f(neq, tn, y, savf)
473 650 yh(i,2) = h*savf(i)
476 IF (nq .EQ. 1) go
to 200
492 700 r = 1.0e0/tesco(2,nqu)
494 710 acor(i) = acor(i)*r
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE * f
subroutine scfode(METH, ELCO, TESCO)
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
charNDArray max(char d, const charNDArray &m)
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
subroutine sstode(NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM, F, JAC, PJAC, SLVS)
charNDArray min(char d, const charNDArray &m)