/**
 * An implementation of ODEX, by E. Hairer and G. Wanner, ported from the Fortran ODEX.F.
 * The original work carries the BSD 2-clause license, and so does this.
 *
 * Copyright (c) 2016 Colin Smith.
 * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following
 * disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the
 * following disclaimer in the documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */ export interface Derivative { ( x , y ) : number [ ] ; } export interface OutputFunction { ( nr , xold , x , y , dense ? ) : boolean | void ; } export enum Outcome { Converged , MaxStepsExceeded , EarlyReturn } export class Solver { n : number ; uRound : number ; maxSteps : number ; initialStepSize : number ; maxStepSize : number ; maxExtrapolationColumns : number ; stepSizeSequence : number ; stabilityCheckCount : number ; stabilityCheckTableLines : number ; denseOutput : boolean ; denseOutputErrorEstimator : boolean ; denseComponents : number [ ] ; interpolationFormulaDegree : number ; stepSizeReductionFactor : number ; stepSizeFac1 : number ; stepSizeFac2 : number ; stepSizeFac3 : number ; stepSizeFac4 : number ; stepSafetyFactor1 : number ; stepSafetyFactor2 : number ; relativeTolerance : number | number [ ] ; absoluteTolerance : number | number [ ] ; debug : boolean ; constructor ( n ) { this . n = n ; this . uRound = 2.3e-16 ; this . maxSteps = 10000 ; this . initialStepSize = 1e-4 ; this . maxStepSize = 0 ; this . maxExtrapolationColumns = 9 ; this . stepSizeSequence = 0 ; this . stabilityCheckCount = 1 ; this . stabilityCheckTableLines = 2 ; this . denseOutput = false ; this . denseOutputErrorEstimator = true ; this . denseComponents = undefined ; this . interpolationFormulaDegree = 4 ; this . stepSizeReductionFactor = 0.5 ; this . stepSizeFac1 = 0.02 ; this . stepSizeFac2 = 4.0 ; this . stepSizeFac3 = 0.8 ; this . stepSizeFac4 = 0.9 ; this . stepSafetyFactor1 = 0.65 ; this . stepSafetyFactor2 = 0.94 ; this . relativeTolerance = 1e-5 ; this . absoluteTolerance = 1e-5 ; this . debug = false ; } grid ( dt , out ) { if ( ! this . denseOutput ) throw new Error ( 'Must set .denseOutput to true when using grid' ) ; let components = this . denseComponents ; if ( ! components ) { components = [ ] ; for ( let i = 0 ; i < this . n ; ++ i ) components . push ( i ) ; } let t ; return ( n , xOld , x , y , interpolate ) => { if ( n === 1 ) { let v = out ( x , y ) ; t = x + dt ; return v ; } while ( t <= x ) { let yf = [ ] ; for ( let i of components ) { yf . push ( interpolate ( i , t ) ) ; } let v = out ( t , yf ) ; if ( v === false ) return false ; t += dt ; } } ; } private static dim = ( n ) => Array ( n + 1 ) ; private static log10 = ( x ) => Math . log ( x ) / Math . LN10 ; private static dim2 ( r , c ) { let a = new Array ( r + 1 ) ; for ( let i = 1 ; i <= r ; ++ i ) a [ i ] = Solver . dim ( c ) ; return a ; } static stepSizeSequence ( nSeq , n ) { const a = new Array ( n + 1 ) ; a [ 0 ] = 0 ; switch ( nSeq ) { case 1 : for ( let i = 1 ; i <= n ; ++ i ) a [ i ] = 2 * i ; break ; case 2 : a [ 1 ] = 2 ; for ( let i = 2 ; i <= n ; ++ i ) a [ i ] = 4 * i - 4 ; break ; case 3 : a [ 1 ] = 2 ; a [ 2 ] = 4 ; a [ 3 ] = 6 ; for ( let i = 4 ; i <= n ; ++ i ) a [ i ] = 2 * a [ i - 2 ] ; break ; case 4 : for ( let i = 1 ; i <= n ; ++ i ) a [ i ] = 4 * i - 2 ; break ; case 5 : for ( let i = 1 ; i <= n ; ++ i ) a [ i ] = 4 * i ; break ; default : throw new Error ( 'invalid stepSizeSequence selected' ) ; } return a ; } solve ( f , x , y0 , xEnd , solOut ? ) { let y = [ 0 ] . concat ( y0 ) ; let dz = Solver . dim ( this . n ) ; let yh1 = Solver . dim ( this . n ) ; let yh2 = Solver . dim ( this . n ) ; if ( this . maxSteps <= 0 ) throw new Error ( 'maxSteps must be positive' ) ; const km = this . maxExtrapolationColumns ; if ( km <= 2 ) throw new Error ( 'maxExtrapolationColumns must be > 2' ) ; const nSeq = this . stepSizeSequence || ( this . denseOutput ? 4 : 1 ) ; if ( nSeq <= 3 && this . denseOutput ) throw new Error ( 'stepSizeSequence incompatible with denseOutput' ) ; if ( this . denseOutput && ! solOut ) throw new Error ( 'denseOutput requires a solution observer function' ) ; if ( this . interpolationFormulaDegree <= 0 || this . interpolationFormulaDegree >= 7 ) throw new Error ( 'bad interpolationFormulaDegree' ) ; let icom = [ 0 ] ; let nrdens = 0 ; if ( this . denseOutput ) { if ( this . denseComponents ) { for ( let c of this . denseComponents ) { if ( c < 0 || c > this . n ) throw new Error ( 'bad dense component: ' + c ) ; icom . push ( c + 1 ) ; ++ nrdens ; } } else { for ( let i = 1 ; i <= this . n ; ++ i ) { icom . push ( i ) ; } nrdens = this . n ; } } if ( this . uRound <= 1e-35 || this . uRound > 1 ) throw new Error ( 'suspicious value of uRound' ) ; const hMax = Math . abs ( this . maxStepSize || xEnd - x ) ; const lfSafe = 2 * km * km + km ; function expandToArray ( x , n ) { const tolArray = [ 0 ] ; if ( Array . isArray ( x ) ) { return tolArray . concat ( x ) ; } else { for ( let i = 0 ; i < n ; ++ i ) tolArray . push ( x ) ; return tolArray ; } } const aTol = expandToArray ( this . absoluteTolerance , this . n ) ; const rTol = expandToArray ( this . relativeTolerance , this . n ) ; let [ nEval , nStep , nAccept , nReject ] = [ 0 , 0 , 0 , 0 ] ; const nrd = Math . max ( 1 , nrdens ) ; const ncom = Math . max ( 1 , ( 2 * km + 5 ) * nrdens ) ; const dens = Solver . dim ( ncom ) ; const fSafe = Solver . dim2 ( lfSafe , nrd ) ; const F = ( x , y , yp ) => { let ret = f ( x , y . slice ( 1 ) ) ; for ( let i = 0 ; i < ret . length ; ++ i ) yp [ i + 1 ] = ret [ i ] ; } ; let odxcor = ( ) : Outcome => { let xOldd ; let hhh ; let kmit ; let acceptStep = ( n ) : boolean => { xOld = x ; x += h ; if ( this . denseOutput ) { kmit = 2 * kc - this . interpolationFormulaDegree + 1 ; for ( let i = 1 ; i <= nrd ; ++ i ) dens [ i ] = y [ icom [ i ] ] ; xOldd = xOld ; hhh = h ; for ( let i = 1 ; i <= nrd ; ++ i ) dens [ nrd + i ] = h * dz [ icom [ i ] ] ; let kln = 2 * nrd ; for ( let i = 1 ; i <= nrd ; ++ i ) dens [ kln + i ] = t [ 1 ] [ icom [ i ] ] ; for ( let j = 2 ; j <= kc ; ++ j ) { let dblenj = nj [ j ] ; for ( let l = j ; l >= 2 ; -- l ) { let factor = ( dblenj / nj [ l - 1 ] ) ** 2 - 1 ; for ( let i = 1 ; i <= nrd ; ++ i ) { ySafe [ l - 1 ] [ i ] = ySafe [ l ] [ i ] + ( ySafe [ l ] [ i ] - ySafe [ l - 1 ] [ i ] ) / factor ; } } } let krn = 4 * nrd ; for ( let i = 1 ; i <= nrd ; ++ i ) dens [ krn + i ] = ySafe [ 1 ] [ i ] ; for ( let i = 1 ; i <= n ; ++ i ) yh1 [ i ] = t [ 1 ] [ i ] ; F ( x , yh1 , yh2 ) ; krn = 3 * nrd ; for ( let i = 1 ; i <= nrd ; ++ i ) dens [ krn + i ] = yh2 [ icom [ i ] ] * h ; for ( let kmi = 1 ; kmi <= kmit ; ++ kmi ) { let kbeg = ( kmi + 1 ) / 2 | 0 ; for ( let kk = kbeg ; kk <= kc ; ++ kk ) { let facnj = ( nj [ kk ] / 2 ) ** ( kmi - 1 ) ; iPt = iPoint [ kk + 1 ] - 2 * kk + kmi ; for ( let i = 1 ; i <= nrd ; ++ i ) { ySafe [ kk ] [ i ] = fSafe [ iPt ] [ i ] * facnj ; } } for ( let j = kbeg + 1 ; j <= kc ; ++ j ) { let dblenj = nj [ j ] ; for ( let l = j ; l >= kbeg + 1 ; -- l ) { let factor = ( dblenj / nj [ l - 1 ] ) ** 2 - 1 ; for ( let i = 1 ; i <= nrd ; ++ i ) { ySafe [ l - 1 ] [ i ] = ySafe [ l ] [ i ] + ( ySafe [ l ] [ i ] - ySafe [ l - 1 ] [ i ] ) / factor ; } } } krn = ( kmi + 4 ) * nrd ; for ( let i = 1 ; i <= nrd ; ++ i ) dens [ krn + i ] = ySafe [ kbeg ] [ i ] * h ; if ( kmi === kmit ) continue ; for ( let kk = ( kmi + 2 ) / 2 | 0 ; kk <= kc ; ++ kk ) { let lbeg = iPoint [ kk + 1 ] ; let lend = iPoint [ kk ] + kmi + 1 ; if ( kmi === 1 && nSeq === 4 ) lend += 2 ; let l ; for ( l = lbeg ; l >= lend ; l -= 2 ) { for ( let i = 1 ; i <= nrd ; ++ i ) { fSafe [ l ] [ i ] -= fSafe [ l - 2 ] [ i ] ; } } if ( kmi === 1 && nSeq === 4 ) { l = lend - 2 ; for ( let i = 1 ; i <= nrd ; ++ i ) fSafe [ l ] [ i ] -= dz [ icom [ i ] ] ; } } for ( let kk = ( kmi + 2 ) / 2 | 0 ; kk <= kc ; ++ kk ) { let lbeg = iPoint [ kk + 1 ] - 1 ; let lend = iPoint [ kk ] + kmi + 2 ; for ( let l = lbeg ; l >= lend ; l -= 2 ) { for ( let i = 1 ; i <= nrd ; ++ i ) { fSafe [ l ] [ i ] -= fSafe [ l - 2 ] [ i ] ; } } } } interp ( nrd , dens , kmit ) ; if ( this . denseOutputErrorEstimator && kmit >= 1 ) { let errint = 0 ; for ( let i = 1 ; i <= nrd ; ++ i ) errint += ( dens [ ( kmit + 4 ) * nrd + i ] / scal [ icom [ i ] ] ) ** 2 ; errint = Math . sqrt ( errint / nrd ) * errfac [ kmit ] ; hoptde = h / Math . max ( errint ** ( 1 / ( kmit + 4 ) ) , 0.01 ) ; if ( errint > 10 ) { h = hoptde ; x = xOld ; ++ nReject ; reject = true ; return true ; } } for ( let i = 1 ; i <= n ; ++ i ) dz [ i ] = yh2 [ i ] ; } for ( let i = 1 ; i <= n ; ++ i ) y [ i ] = t [ 1 ] [ i ] ; ++ nAccept ; if ( solOut ) { if ( solOut ( nAccept + 1 , xOld , x , y . slice ( 1 ) , this . denseOutput && contex ( xOldd , hhh , kmit , dens , icom ) ) === false ) return false ; } let kopt ; if ( kc === 2 ) { kopt = Math . min ( 3 , km - 1 ) ; if ( reject ) kopt = 2 ; } else { if ( kc <= k ) { kopt = kc ; if ( w [ kc - 1 ] < w [ kc ] * this . stepSizeFac3 ) kopt = kc - 1 ; if ( w [ kc ] < w [ kc - 1 ] * this . stepSizeFac4 ) kopt = Math . min ( kc + 1 , km - 1 ) ; } else { kopt = kc - 1 ; if ( kc > 3 && w [ kc - 2 ] < w [ kc - 1 ] * this . stepSizeFac3 ) kopt = kc - 2 ; if ( w [ kc ] < w [ kopt ] * this . stepSizeFac4 ) kopt = Math . min ( kc , km - 1 ) ; } } if ( reject ) { k = Math . min ( kopt , kc ) ; h = posneg * Math . min ( Math . abs ( h ) , Math . abs ( hh [ k ] ) ) ; reject = false ; return true ; } if ( kopt <= kc ) { h = hh [ kopt ] ; } else { if ( kc < k && w [ kc ] < w [ kc - 1 ] * this . stepSizeFac4 ) { h = hh [ kc ] * a [ kopt + 1 ] / a [ kc ] ; } else { h = hh [ kc ] * a [ kopt ] / a [ kc ] ; } } k = kopt ; h = posneg * Math . abs ( h ) ; return true ; } ; let midex = ( j ) : void => { const dy = Solver . dim ( this . n ) ; const hj = h / nj [ j ] ; for ( let i = 1 ; i <= this . n ; ++ i ) { yh1 [ i ] = y [ i ] ; yh2 [ i ] = y [ i ] + hj * dz [ i ] ; } const m = nj [ j ] - 1 ; const njMid = ( nj [ j ] / 2 ) | 0 ; for ( let mm = 1 ; mm <= m ; ++ mm ) { if ( this . denseOutput && mm === njMid ) { for ( let i = 1 ; i <= nrd ; ++ i ) { ySafe [ j ] [ i ] = yh2 [ icom [ i ] ] ; } } F ( x + hj * mm , yh2 , dy ) ; if ( this . denseOutput && Math . abs ( mm - njMid ) <= 2 * j - 1 ) { ++ iPt ; for ( let i = 1 ; i <= nrd ; ++ i ) { fSafe [ iPt ] [ i ] = dy [ icom [ i ] ] ; } } for ( let i = 1 ; i <= this . n ; ++ i ) { let ys = yh1 [ i ] ; yh1 [ i ] = yh2 [ i ] ; yh2 [ i ] = ys + 2 * hj * dy [ i ] ; } if ( mm <= this . stabilityCheckCount && j <= this . stabilityCheckTableLines ) { let del1 = 0 ; for ( let i = 1 ; i <= this . n ; ++ i ) { del1 += ( dz [ i ] / scal [ i ] ) ** 2 ; } let del2 = 0 ; for ( let i = 1 ; i <= this . n ; ++ i ) { del2 += ( ( dy [ i ] - dz [ i ] ) / scal [ i ] ) ** 2 ; } const quot = del2 / Math . max ( this . uRound , del1 ) ; if ( quot > 4 ) { ++ nEval ; atov = true ; h *= this . stepSizeReductionFactor ; reject = true ; return ; } } } F ( x + h , yh2 , dy ) ; if ( this . denseOutput && njMid <= 2 * j - 1 ) { ++ iPt ; for ( let i = 1 ; i <= nrd ; ++ i ) { fSafe [ iPt ] [ i ] = dy [ icom [ i ] ] ; } } for ( let i = 1 ; i <= this . n ; ++ i ) { t [ j ] [ i ] = ( yh1 [ i ] + yh2 [ i ] + hj * dy [ i ] ) / 2 ; } nEval += nj [ j ] ; if ( j === 1 ) return ; const dblenj = nj [ j ] ; let fac ; for ( let l = j ; l > 1 ; -- l ) { fac = ( dblenj / nj [ l - 1 ] ) ** 2 - 1 ; for ( let i = 1 ; i <= this . n ; ++ i ) { t [ l - 1 ] [ i ] = t [ l ] [ i ] + ( t [ l ] [ i ] - t [ l - 1 ] [ i ] ) / fac ; } } err = 0 ; for ( let i = 1 ; i <= this . n ; ++ i ) { let t1i = Math . max ( Math . abs ( y [ i ] ) , Math . abs ( t [ 1 ] [ i ] ) ) ; scal [ i ] = aTol [ i ] + rTol [ i ] * t1i ; err += ( ( t [ 1 ] [ i ] - t [ 2 ] [ i ] ) / scal [ i ] ) ** 2 ; } err = Math . sqrt ( err / this . n ) ; if ( err * this . uRound >= 1 || ( j > 2 && err >= errOld ) ) { atov = true ; h *= this . stepSizeReductionFactor ; reject = true ; return ; } errOld = Math . max ( 4 * err , 1 ) ; let exp0 = 1 / ( 2 * j - 1 ) ; let facMin = this . stepSizeFac1 ** exp0 ; fac = Math . min ( this . stepSizeFac2 / facMin , Math . max ( facMin , ( err / this . stepSafetyFactor1 ) ** exp0 / this . stepSafetyFactor2 ) ) ; fac = 1 / fac ; hh [ j ] = Math . min ( Math . abs ( h ) * fac , hMax ) ; w [ j ] = a [ j ] / hh [ j ] ; } ; const interp = ( n , y , imit ) => { let a = new Array ( 31 ) ; for ( let i = 1 ; i <= n ; ++ i ) { let y0 = y [ i ] ; let y1 = y [ 2 * n + i ] ; let yp0 = y [ n + i ] ; let yp1 = y [ 3 * n + i ] ; let yDiff = y1 - y0 ; let aspl = - yp1 + yDiff ; let bspl = yp0 - yDiff ; y [ n + i ] = yDiff ; y [ 2 * n + i ] = aspl ; y [ 3 * n + i ] = bspl ; if ( imit < 0 ) continue ; let ph0 = ( y0 + y1 ) * 0.5 + 0.125 * ( aspl + bspl ) ; let ph1 = yDiff + ( aspl - bspl ) * 0.25 ; let ph2 = - ( yp0 - yp1 ) ; let ph3 = 6 * ( bspl - aspl ) ; if ( imit >= 1 ) { a [ 1 ] = 16 * ( y [ 5 * n + i ] - ph1 ) ; if ( imit >= 3 ) { a [ 3 ] = 16 * ( y [ 7 * n + i ] - ph3 + 3 * a [ 1 ] ) ; if ( imit >= 5 ) { for ( let im = 5 ; im <= imit ; im += 2 ) { let fac1 = im * ( im - 1 ) / 2 ; let fac2 = fac1 * ( im - 2 ) * ( im - 3 ) * 2 ; a [ im ] = 16 * ( y [ ( im + 4 ) * n + i ] + fac1 * a [ im - 2 ] - fac2 * a [ im - 4 ] ) ; } } } } a [ 0 ] = ( y [ 4 * n + i ] - ph0 ) * 16 ; if ( imit >= 2 ) { a [ 2 ] = ( y [ n * 6 + i ] - ph2 + a [ 0 ] ) * 16 ; if ( imit >= 4 ) { for ( let im = 4 ; im <= imit ; im += 2 ) { let fac1 = im * ( im - 1 ) / 2 ; let fac2 = im * ( im - 1 ) * ( im - 2 ) * ( im - 3 ) ; a [ im ] = ( y [ n * ( im + 4 ) + i ] + a [ im - 2 ] * fac1 - a [ im - 4 ] * fac2 ) * 16 ; } } } for ( let im = 0 ; im <= imit ; ++ im ) y [ n * ( im + 4 ) + i ] = a [ im ] ; } } ; const contex = ( xOld , h , imit , y , icom ) => { return ( c , x ) => { let i = 0 ; for ( let j = 1 ; j <= nrd ; ++ j ) { if ( icom [ j ] === c + 1 ) i = j ; } if ( i === 0 ) throw new Error ( 'no dense output available for component ' + c ) ; const theta = ( x - xOld ) / h ; const theta1 = 1 - theta ; const phthet = y [ i ] + theta * ( y [ nrd + i ] + theta1 * ( y [ 2 * nrd + i ] * theta + y [ 3 * nrd + i ] * theta1 ) ) ; if ( imit < 0 ) return phthet ; const thetah = theta - 0.5 ; let ret = y [ nrd * ( imit + 4 ) + i ] ; for ( let im = imit ; im >= 1 ; -- im ) { ret = y [ nrd * ( im + 3 ) + i ] + ret * thetah / im ; } return phthet + ( theta * theta1 ) ** 2 * ret ; } ; } ; const ySafe = Solver . dim2 ( km , nrd ) ; const hh = Solver . dim ( km ) ; const t = Solver . dim2 ( km , this . n ) ; const nj = Solver . stepSizeSequence ( nSeq , km ) ; const a = Solver . dim ( km ) ; a [ 1 ] = 1 + nj [ 1 ] ; for ( let i = 2 ; i <= km ; ++ i ) { a [ i ] = a [ i - 1 ] + nj [ i ] ; } const scal = Solver . dim ( this . n ) ; for ( let i = 1 ; i <= this . n ; ++ i ) { scal [ i ] = aTol [ i ] + rTol [ i ] + Math . abs ( y [ i ] ) ; } const posneg = xEnd - x >= 0 ? 1 : - 1 ; let k = Math . max ( 2 , Math . min ( km - 1 , Math . floor ( - Solver . log10 ( rTol [ 1 ] + 1e-40 ) * 0.6 + 1.5 ) ) ) ; let h = Math . max ( Math . abs ( this . initialStepSize ) , 1e-4 ) ; h = posneg * Math . min ( h , hMax , Math . abs ( xEnd - x ) / 2 ) ; const iPoint = Solver . dim ( km + 1 ) ; const errfac = Solver . dim ( 2 * km ) ; let xOld = x ; let iPt = 0 ; if ( solOut ) { if ( this . denseOutput ) { iPoint [ 1 ] = 0 ; for ( let i = 1 ; i <= km ; ++ i ) { let njAdd = 4 * i - 2 ; if ( nj [ i ] > njAdd ) ++ njAdd ; iPoint [ i + 1 ] = iPoint [ i ] + njAdd ; } for ( let mu = 1 ; mu <= 2 * km ; ++ mu ) { let errx = Math . sqrt ( mu / ( mu + 4 ) ) * 0.5 ; let prod = ( 1 / ( mu + 4 ) ) ** 2 ; for ( let j = 1 ; j <= mu ; ++ j ) prod *= errx / j ; errfac [ mu ] = prod ; } iPt = 0 ; } if ( false === solOut ( nAccept + 1 , xOld , x , y . slice ( 1 ) ) ) { return Outcome . EarlyReturn ; } } let err = 0 ; let errOld = 1e10 ; let hoptde = posneg * hMax ; const w = Solver . dim ( km ) ; w [ 1 ] = 0 ; let reject = false ; let last = false ; let atov ; let kc = 0 ; enum STATE { Start , BasicIntegrationStep , ConvergenceStep , HopeForConvergence , Accept , Reject } let state = STATE . Start ; loop : while ( true ) { this . debug && console . log ( 'STATE' , STATE [ state ] , nStep , xOld , x , h , k , kc , hoptde ) ; switch ( state ) { case STATE . Start : atov = false ; if ( 0.1 * Math . abs ( xEnd - x ) <= Math . abs ( x ) * this . uRound ) break loop ; h = posneg * Math . min ( Math . abs ( h ) , Math . abs ( xEnd - x ) , hMax , Math . abs ( hoptde ) ) ; if ( ( x + 1.01 * h - xEnd ) * posneg > 0 ) { h = xEnd - x ; last = true ; } if ( nStep === 0 || ! this . denseOutput ) { F ( x , y , dz ) ; ++ nEval ; } if ( nStep === 0 || last ) { iPt = 0 ; ++ nStep ; for ( let j = 1 ; j <= k ; ++ j ) { kc = j ; midex ( j ) ; if ( atov ) continue loop ; if ( j > 1 && err <= 1 ) { state = STATE . Accept ; continue loop ; } } state = STATE . HopeForConvergence ; continue ; } state = STATE . BasicIntegrationStep ; continue ; case STATE . BasicIntegrationStep : iPt = 0 ; ++ nStep ; if ( nStep >= this . maxSteps ) { return Outcome . MaxStepsExceeded ; } kc = k - 1 ; for ( let j = 1 ; j <= kc ; ++ j ) { midex ( j ) ; if ( atov ) { state = STATE . Start ; continue loop ; } } if ( k === 2 || reject ) { state = STATE . ConvergenceStep ; } else { if ( err <= 1 ) { state = STATE . Accept ; } else if ( err > ( ( nj [ k + 1 ] * nj [ k ] ) / 4 ) ** 2 ) { state = STATE . Reject ; } else state = STATE . ConvergenceStep ; } continue ; case STATE . ConvergenceStep : midex ( k ) ; if ( atov ) { state = STATE . Start ; continue ; } kc = k ; if ( err <= 1 ) { state = STATE . Accept ; continue ; } state = STATE . HopeForConvergence ; continue ; case STATE . HopeForConvergence : if ( err > ( nj [ k + 1 ] / 2 ) ** 2 ) { state = STATE . Reject ; continue ; } kc = k + 1 ; midex ( kc ) ; if ( atov ) state = STATE . Start ; else if ( err > 1 ) state = STATE . Reject ; else state = STATE . Accept ; continue ; case STATE . Accept : if ( ! acceptStep ( this . n ) ) return Outcome . EarlyReturn ; state = STATE . Start ; continue ; case STATE . Reject : k = Math . min ( k , kc , km - 1 ) ; if ( k > 2 && w [ k - 1 ] < w [ k ] * this . stepSizeFac3 ) k -= 1 ; ++ nReject ; h = posneg * hh [ k ] ; reject = true ; state = STATE . BasicIntegrationStep ; } } return Outcome . Converged ; } ; const outcome = odxcor ( ) ; return { y : y . slice ( 1 ) , outcome : outcome , nStep : nStep , xEnd : xEnd , nAccept : nAccept , nReject : nReject , nEval : nEval } ; } }