//+------------------------------------------------------------------+
//|                                                diffequations.mqh |
//|            Copyright 2003-2022 Sergey Bochkanov (ALGLIB project) |
//|                             Copyright 2012-2023, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
//| Implementation of ALGLIB library in MetaQuotes Language 5        |
//|                                                                  |
//| The features of the library include:                             |
//| - Linear algebra (direct algorithms, EVD, SVD)                   |
//| - Solving systems of linear and non-linear equations             |
//| - Interpolation                                                  |
//| - Optimization                                                   |
//| - FFT (Fast Fourier Transform)                                   |
//| - Numerical integration                                          |
//| - Linear and nonlinear least-squares fitting                     |
//| - Ordinary differential equations                                |
//| - Computation of special functions                               |
//| - Descriptive statistics and hypothesis testing                  |
//| - Data analysis - classification, regression                     |
//| - Implementing linear algebra algorithms, interpolation, etc.    |
//|   in high-precision arithmetic (using MPFR)                      |
//|                                                                  |
//| This file is free software; you can redistribute it and/or       |
//| modify it under the terms of the GNU General Public License as   |
//| published by the Free Software Foundation (www.fsf.org); either  |
//| version 2 of the License, or (at your option) any later version. |
//|                                                                  |
//| This program is distributed in the hope that it will be useful,  |
//| but WITHOUT ANY WARRANTY; without even the implied warranty of   |
//| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the     |
//| GNU General Public License for more details.                     |
//+------------------------------------------------------------------+
#include "matrix.mqh"
#include "alglibinternal.mqh"
//+------------------------------------------------------------------+
//| Auxiliary class for CODESolver                                   |
//+------------------------------------------------------------------+
class CODESolverState
  {
public:
   int               m_n;
   int               m_m;
   double            m_xscale;
   double            m_h;
   double            m_eps;
   bool              m_fraceps;
   int               m_repterminationtype;
   int               m_repnfev;
   int               m_solvertype;
   bool              m_needdy;
   double            m_x;
   RCommState        m_rstate;
   //--- arrays
   double            m_yc[];
   double            m_escale[];
   double            m_xg[];
   double            m_y[];
   double            m_dy[];
   double            m_yn[];
   double            m_yns[];
   double            m_rka[];
   double            m_rkc[];
   double            m_rkcs[];
   //--- matrices
   CMatrixDouble     m_ytbl;
   CMatrixDouble     m_rkb;
   CMatrixDouble     m_rkk;

public:
                     CODESolverState(void);
                    ~CODESolverState(void) {}

   void              Copy(CODESolverState &obj);
  };
//+------------------------------------------------------------------+
//| Constructor                                                      |
//+------------------------------------------------------------------+
CODESolverState::CODESolverState(void)
  {
   m_n=0;
   m_m=0;
   m_xscale=0;
   m_h=0;
   m_eps=0;
   m_fraceps=false;
   m_repterminationtype=0;
   m_repnfev=0;
   m_solvertype=0;
   m_needdy=false;
   m_x=0;
  }
//+------------------------------------------------------------------+
//| Copy                                                             |
//+------------------------------------------------------------------+
void CODESolverState::Copy(CODESolverState &obj)
  {
//--- copy variables
   m_n=obj.m_n;
   m_m=obj.m_m;
   m_xscale=obj.m_xscale;
   m_h=obj.m_h;
   m_eps=obj.m_eps;
   m_fraceps=obj.m_fraceps;
   m_repterminationtype=obj.m_repterminationtype;
   m_repnfev=obj.m_repnfev;
   m_solvertype=obj.m_solvertype;
   m_needdy=obj.m_needdy;
   m_x=obj.m_x;
   m_rstate.Copy(obj.m_rstate);
//--- copy arrays
   ArrayCopy(m_yc,obj.m_yc);
   ArrayCopy(m_escale,obj.m_escale);
   ArrayCopy(m_xg,obj.m_xg);
   ArrayCopy(m_y,obj.m_y);
   ArrayCopy(m_dy,obj.m_dy);
   ArrayCopy(m_yn,obj.m_yn);
   ArrayCopy(m_yns,obj.m_yns);
   ArrayCopy(m_rka,obj.m_rka);
   ArrayCopy(m_rkc,obj.m_rkc);
   ArrayCopy(m_rkcs,obj.m_rkcs);
//--- copy matrices
   m_ytbl=obj.m_ytbl;
   m_rkb=obj.m_rkb;
   m_rkk=obj.m_rkk;
  }
//+------------------------------------------------------------------+
//| This class is a shell for class CODESolverState                  |
//+------------------------------------------------------------------+
class CODESolverStateShell
  {
private:
   CODESolverState   m_innerobj;

public:
   //--- constructors, destructor
                     CODESolverStateShell(void) {}
                     CODESolverStateShell(CODESolverState &obj) { m_innerobj.Copy(obj); }
                    ~CODESolverStateShell(void) {}
   //--- methods
   bool              GetNeedDY(void);
   void              SetNeedDY(const bool b);
   double            GetX(void);
   void              SetX(const double d);
   CODESolverState  *GetInnerObj(void);
  };
//+------------------------------------------------------------------+
//| Returns the value of the variable needdy                         |
//+------------------------------------------------------------------+
bool CODESolverStateShell::GetNeedDY(void)
  {
   return(m_innerobj.m_needdy);
  }
//+------------------------------------------------------------------+
//| Changing the value of the variable needdy                        |
//+------------------------------------------------------------------+
void CODESolverStateShell::SetNeedDY(const bool b)
  {
   m_innerobj.m_needdy=b;
  }
//+------------------------------------------------------------------+
//| Returns the value of the variable x                              |
//+------------------------------------------------------------------+
double CODESolverStateShell::GetX(void)
  {
   return(m_innerobj.m_x);
  }
//+------------------------------------------------------------------+
//| Changing the value of the variable x                             |
//+------------------------------------------------------------------+
void CODESolverStateShell::SetX(const double d)
  {
   m_innerobj.m_x=d;
  }
//+------------------------------------------------------------------+
//| Return object of class                                           |
//+------------------------------------------------------------------+
CODESolverState *CODESolverStateShell::GetInnerObj(void)
  {
   return(GetPointer(m_innerobj));
  }
//+------------------------------------------------------------------+
//| Auxiliary class for CODESolver                                   |
//+------------------------------------------------------------------+
class CODESolverReport
  {
public:
   //--- class variables
   int               m_nfev;
   int               m_terminationtype;
   //--- constructor, destructor
                     CODESolverReport(void) { ZeroMemory(this); }
                    ~CODESolverReport(void) {}
   //--- copy
   void              Copy(CODESolverReport &obj);
  };
//+------------------------------------------------------------------+
//| Copy                                                             |
//+------------------------------------------------------------------+
void CODESolverReport::Copy(CODESolverReport &obj)
  {
//--- copy variables
   m_nfev=obj.m_nfev;
   m_terminationtype=obj.m_terminationtype;
  }
//+------------------------------------------------------------------+
//| This class is a shell for class CODESolverReport                 |
//+------------------------------------------------------------------+
class CODESolverReportShell
  {
private:
   CODESolverReport  m_innerobj;

public:
   //--- constructor, destructor
                     CODESolverReportShell(void) {}
                     CODESolverReportShell(CODESolverReport &obj) { m_innerobj.Copy(obj); }
                    ~CODESolverReportShell(void) {}
   //--- methods
   int               GetNFev(void);
   void              SetNFev(const int i);
   int               GetTerminationType(void);
   void              SetTerminationType(const int i);
   CODESolverReport *GetInnerObj(void);
  };
//+------------------------------------------------------------------+
//| Returns the value of the variable nfev                           |
//+------------------------------------------------------------------+
int CODESolverReportShell::GetNFev(void)
  {
   return(m_innerobj.m_nfev);
  }
//+------------------------------------------------------------------+
//| Changing the value of the variable nfev                          |
//+------------------------------------------------------------------+
void CODESolverReportShell::SetNFev(const int i)
  {
   m_innerobj.m_nfev=i;
  }
//+------------------------------------------------------------------+
//| Returns the value of the variable terminationtype                |
//+------------------------------------------------------------------+
int CODESolverReportShell::GetTerminationType(void)
  {
   return(m_innerobj.m_terminationtype);
  }
//+------------------------------------------------------------------+
//| Changing the value of the variable terminationtype               |
//+------------------------------------------------------------------+
void CODESolverReportShell::SetTerminationType(const int i)
  {
   m_innerobj.m_terminationtype=i;
  }
//+------------------------------------------------------------------+
//| Return object of class                                           |
//+------------------------------------------------------------------+
CODESolverReport *CODESolverReportShell::GetInnerObj(void)
  {
   return(GetPointer(m_innerobj));
  }
//+------------------------------------------------------------------+
//| Solution of ordinary differential equations                      |
//+------------------------------------------------------------------+
class CODESolver
  {
public:
   //--- class constants
   static const double m_odesolvermaxgrow;
   static const double m_odesolvermaxshrink;
   //--- public methods
   static void       ODESolverRKCK(double &y[],const int n,double &x[],const int m,const double eps,const double h,CODESolverState &state);
   static void       ODESolverResults(CODESolverState &state,int &m,double &xtbl[],CMatrixDouble &ytbl,CODESolverReport &rep);
   static bool       ODESolverIteration(CODESolverState &state);

private:
   //--- private method
   static void       ODESolverInit(int solvertype,double &y[],const int n,double &x[],const int m,const double eps,double h,CODESolverState &state);
   //--- auxiliary functions for ODESolverIteration
   static void       Func_lbl_rcomm(CODESolverState &state,int n,int m,int i,int j,int k,int klimit,bool gridpoint,double xc,double v,double h,double h2,double err,double maxgrowpow);
   static bool       Func_lbl_6(CODESolverState &state,int &n,int &m,int &i,int &j,int &k,int &klimit,bool &gridpoint,double &xc,double &v,double &h,double &h2,double &err,double &maxgrowpow);
   static bool       Func_lbl_8(CODESolverState &state,int &n,int &m,int &i,int &j,int &k,int &klimit,bool &gridpoint,double &xc,double &v,double &h,double &h2,double &err,double &maxgrowpow);
   static bool       Func_lbl_10(CODESolverState &state,int &n,int &m,int &i,int &j,int &k,int &klimit,bool &gridpoint,double &xc,double &v,double &h,double &h2,double &err,double &maxgrowpow);
  };
//+------------------------------------------------------------------+
//| Initialize constants                                             |
//+------------------------------------------------------------------+
const double CODESolver::m_odesolvermaxgrow=3.0;
const double CODESolver::m_odesolvermaxshrink=10.0;
//+------------------------------------------------------------------+
//| Cash-Karp adaptive ODE solver.                                   |
//| This subroutine solves ODE  Y'=f(Y,x) with initial conditions    |
//| Y(xs)=Ys (here Y may be single variable or vector of N variables)|
//| INPUT PARAMETERS:                                                |
//|     Y       -   initial conditions, array[0..N-1].               |
//|                 contains values of Y[] at X[0]                   |
//|     N       -   system size                                      |
//|     X       -   points at which Y should be tabulated,           |
//|                 array[0..M-1] integrations starts at X[0], ends  |
//|                 at X[M-1], intermediate values at X[i] are       |
//|                 returned too.                                    |
//|                 SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!!|
//|     M       -   number of intermediate points + first point +    |
//|                 last point:                                      |
//|                 * M>2 means that you need both Y(X[M-1]) and M-2 |
//|                   values at intermediate points                  |
//|                 * M=2 means that you want just to integrate from |
//|                   X[0] to X[1] and don't interested in           |
//|                   intermediate values.                           |
//|                 * M=1 means that you don't want to integrate :)  |
//|                   it is degenerate case, but it will be handled  |
//|                   correctly.                                     |
//|                 * M<1 means error                                |
//|     Eps     -   tolerance (absolute/relative error on each step  |
//|                 will be less than Eps). When passing:            |
//|                 * Eps>0, it means desired ABSOLUTE error         |
//|                 * Eps<0, it means desired RELATIVE error.        |
//|                   Relative errors are calculated with respect to |
//|                   maximum values of Y seen so far. Be careful to |
//|                   use this criterion when starting from Y[] that |
//|                   are close to zero.                             |
//|     H       -   initial step lenth, it will be adjusted          |
//|                 automatically after the first step. If H=0, step |
//|                 will be selected automatically (usualy it will   |
//|                 be equal to 0.001 of min(x[i]-x[j])).            |
//| OUTPUT PARAMETERS                                                |
//|     State   -   structure which stores algorithm state between   |
//|                 subsequent calls of OdeSolverIteration. Used     |
//|                 for reverse communication. This structure should |
//|                 be passed to the OdeSolverIteration subroutine.  |
//| SEE ALSO                                                         |
//|     AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults|
//+------------------------------------------------------------------+
void CODESolver::ODESolverRKCK(double &y[],const int n,double &x[],
                               const int m,const double eps,
                               const double h,CODESolverState &state)
  {
//--- check
   if(!CAp::Assert(n>=1,"ODESolverRKCK: N<1!"))
      return;
//--- check
   if(!CAp::Assert(m>=1,"ODESolverRKCK: M<1!"))
      return;
//--- check
   if(!CAp::Assert(CAp::Len(y)>=n,"ODESolverRKCK: Length(Y)<N!"))
      return;
//--- check
   if(!CAp::Assert(CAp::Len(x)>=m,"ODESolverRKCK: Length(X)<M!"))
      return;
//--- check
   if(!CAp::Assert(CApServ::IsFiniteVector(y,n),"ODESolverRKCK: Y contains infinite or NaN values!"))
      return;
//--- check
   if(!CAp::Assert(CApServ::IsFiniteVector(x,m),"ODESolverRKCK: Y contains infinite or NaN values!"))
      return;
//--- check
   if(!CAp::Assert(CMath::IsFinite(eps),"ODESolverRKCK: Eps is not finite!"))
      return;
//--- check
   if(!CAp::Assert(eps!=0.0,"ODESolverRKCK: Eps is zero!"))
      return;
//--- check
   if(!CAp::Assert(CMath::IsFinite(h),"ODESolverRKCK: H is not finite!"))
      return;
//--- function call
   ODESolverInit(0,y,n,x,m,eps,h,state);
  }
//+------------------------------------------------------------------+
//| ODE solver results                                               |
//| Called after OdeSolverIteration returned False.                  |
//| INPUT PARAMETERS:                                                |
//|     State   -   algorithm state (used by OdeSolverIteration).    |
//| OUTPUT PARAMETERS:                                               |
//|     M       -   number of tabulated values, M>=1                 |
//|     XTbl    -   array[0..M-1], values of X                       |
//|     YTbl    -   array[0..M-1,0..N-1], values of Y in X[i]        |
//|     Rep     -   solver report:                                   |
//|                 * Rep.TerminationType completetion code:         |
//|                     * -2    X is not ordered  by                 |
//|                             ascending/descending or there are    |
//|                             non-distinct X[], i.e.  X[i]=X[i+1]  |
//|                     * -1    incorrect parameters were specified  |
//|                     *  1    task has been solved                 |
//|                 * Rep.NFEV contains number of function           |
//|                   calculations                                   |
//+------------------------------------------------------------------+
void CODESolver::ODESolverResults(CODESolverState &state,int &m,
                                  double &xtbl[],CMatrixDouble &ytbl,
                                  CODESolverReport &rep)
  {
//--- create variables
   double v=0;
   int    i=0;
   int    i_=0;
//--- initialization
   m=0;
   rep.m_terminationtype=state.m_repterminationtype;
//--- check
   if(rep.m_terminationtype>0)
     {
      //--- change values
      m=state.m_m;
      rep.m_nfev=state.m_repnfev;
      //--- allocation
      ArrayResize(xtbl,state.m_m);
      v=state.m_xscale;
      //--- calculation
      for(i_=0; i_<=state.m_m-1; i_++)
         xtbl[i_]=v*state.m_xg[i_];
      //--- allocation
      ytbl.Resize(state.m_m,state.m_n);
      for(i=0; i<=state.m_m-1; i++)
        {
         for(i_=0; i_<=state.m_n-1; i_++)
            ytbl.Set(i,i_,state.m_ytbl[i][i_]);
        }
     }
   else
      rep.m_nfev=0;
  }
//+------------------------------------------------------------------+
//| Internal initialization subroutine                               |
//+------------------------------------------------------------------+
void CODESolver::ODESolverInit(int solvertype,double &y[],const int n,
                               double &x[],const int m,const double eps,
                               double h,CODESolverState &state)
  {
//--- create variables
   int    i=0;
   double v=0;
   int    i_=0;
//--- Prepare RComm
   state.m_rstate.ia.Resize(6);
   ArrayResizeAL(state.m_rstate.ba,1);
   state.m_rstate.ra.Resize(6);
   state.m_rstate.stage=-1;
   state.m_needdy=false;
//--- check parameters.
   if((n<=0 || m<1) || eps==0.0)
     {
      state.m_repterminationtype=-1;
      return;
     }
//--- check
   if(h<0.0)
      h=-h;
//--- quick exit if necessary.
//--- after this block we assume that M>1
   if(m==1)
     {
      //--- change values
      state.m_repnfev=0;
      state.m_repterminationtype=1;
      state.m_ytbl.Resize(1,n);
      for(i_=0; i_<=n-1; i_++)
         state.m_ytbl.Set(0,i_,y[i_]);
      //--- allocation
      ArrayResize(state.m_xg,m);
      for(i_=0; i_<=m-1; i_++)
         state.m_xg[i_]=x[i_];
      //--- exit the function
      return;
     }
//--- check again: correct order of X[]
   if(x[1]==x[0])
     {
      state.m_repterminationtype=-2;
      return;
     }
   for(i=1; i<=m-1; i++)
     {
      //--- check
      if((x[1]>x[0] && x[i]<=x[i-1]) || (x[1]<x[0] && x[i]>=x[i-1]))
        {
         state.m_repterminationtype=-2;
         return;
        }
     }
//--- auto-select H if necessary
   if(h==0.0)
     {
      v=MathAbs(x[1]-x[0]);
      for(i=2; i<=m-1; i++)
         v=MathMin(v,MathAbs(x[i]-x[i-1]));
      h=0.001*v;
     }
//--- store parameters
   state.m_n=n;
   state.m_m=m;
   state.m_h=h;
   state.m_eps=MathAbs(eps);
   state.m_fraceps=eps<0.0;
//--- allocation
   ArrayResize(state.m_xg,m);
   for(i_=0; i_<=m-1; i_++)
      state.m_xg[i_]=x[i_];
//--- check
   if(x[1]>x[0])
      state.m_xscale=1;
   else
     {
      state.m_xscale=-1;
      for(i_=0; i_<=m-1; i_++)
         state.m_xg[i_]=-1*state.m_xg[i_];
     }
//--- allocation
   ArrayResize(state.m_yc,n);
   for(i_=0; i_<=n-1; i_++)
      state.m_yc[i_]=y[i_];
//--- change values
   state.m_solvertype=solvertype;
   state.m_repterminationtype=0;
//--- Allocate arrays
   ArrayResize(state.m_y,n);
   ArrayResize(state.m_dy,n);
  }
//+------------------------------------------------------------------+
//| Iterative method                                                 |
//+------------------------------------------------------------------+
bool CODESolver::ODESolverIteration(CODESolverState &state)
  {
//--- create variables
   int    n=0;
   int    m=0;
   int    i=0;
   int    j=0;
   int    k=0;
   double xc=0;
   double v=0;
   double h=0;
   double h2=0;
   bool   gridpoint;
   double err=0;
   double maxgrowpow=0;
   int    klimit=0;
   int    i_=0;
//--- This code initializes locals by:
//--- * random values determined during code
//---   generation - on first subroutine call
//--- * values from previous call - on subsequent calls
   if(state.m_rstate.stage>=0)
     {
      //--- initialization
      n=state.m_rstate.ia[0];
      m=state.m_rstate.ia[1];
      i=state.m_rstate.ia[2];
      j=state.m_rstate.ia[3];
      k=state.m_rstate.ia[4];
      klimit=state.m_rstate.ia[5];
      gridpoint=state.m_rstate.ba[0];
      xc=state.m_rstate.ra[0];
      v=state.m_rstate.ra[1];
      h=state.m_rstate.ra[2];
      h2=state.m_rstate.ra[3];
      err=state.m_rstate.ra[4];
      maxgrowpow=state.m_rstate.ra[5];
     }
   else
     {
      //--- initialization
      n=-983;
      m=-989;
      i=-834;
      j=900;
      k=-287;
      klimit=364;
      gridpoint=false;
      xc=-338;
      v=-686;
      h=912;
      h2=585;
      err=497;
      maxgrowpow=-271;
     }
//--- check
   if(state.m_rstate.stage==0)
     {
      //--- change values
      state.m_needdy=false;
      state.m_repnfev=state.m_repnfev+1;
      v=h*state.m_xscale;
      for(i_=0; i_<=n-1; i_++)
         state.m_rkk.Set(k,i_,v*state.m_dy[i_]);
      //--- update YN/YNS
      v=state.m_rkc[k];
      for(i_=0; i_<=n-1; i_++)
         state.m_yn[i_]=state.m_yn[i_]+v*state.m_rkk[k][i_];
      v=state.m_rkcs[k];
      for(i_=0; i_<=n-1; i_++)
         state.m_yns[i_]=state.m_yns[i_]+v*state.m_rkk[k][i_];
      k=k+1;
      return(Func_lbl_8(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
     }
//--- Routine body
//--- prepare
   if(state.m_repterminationtype!=0)
      return(false);
//--- change values
   n=state.m_n;
   m=state.m_m;
   h=state.m_h;
   maxgrowpow=MathPow(m_odesolvermaxgrow,5);
   state.m_repnfev=0;
//--- some preliminary checks for internal errors
//--- after this we assume that H>0 and M>1
   if(!CAp::Assert(state.m_h>0.0,"ODESolver: internal error"))
      return(false);
//--- check
   if(!CAp::Assert(m>1,"ODESolverIteration: internal error"))
      return(false);
//--- choose solver
   if(state.m_solvertype!=0)
      return(false);
//--- Cask-Karp solver
//--- Prepare coefficients table.
//--- Check it for errors
   ArrayResize(state.m_rka,6);
//--- calculation
   state.m_rka[0]=0;
   state.m_rka[1]=1.0/5.0;
   state.m_rka[2]=3.0/10.0;
   state.m_rka[3]=3.0/5.0;
   state.m_rka[4]=1;
   state.m_rka[5]=7.0/8.0;
   state.m_rkb.Resize(6,5);
   state.m_rkb.Set(1,0,1.0/5.0);
   state.m_rkb.Set(2,0,3.0/40.0);
   state.m_rkb.Set(2,1,9.0/40.0);
   state.m_rkb.Set(3,0,3.0/10.0);
   state.m_rkb.Set(3,1,-(9.0/10.0));
   state.m_rkb.Set(3,2,6.0/5.0);
   state.m_rkb.Set(4,0,-(11.0/54.0));
   state.m_rkb.Set(4,1,5.0/2.0);
   state.m_rkb.Set(4,2,-(70.0/27.0));
   state.m_rkb.Set(4,3,35.0/27.0);
   state.m_rkb.Set(5,0,1631.0/55296.0);
   state.m_rkb.Set(5,1,175.0/512.0);
   state.m_rkb.Set(5,2,575.0/13824.0);
   state.m_rkb.Set(5,3,44275.0/110592.0);
   state.m_rkb.Set(5,4,253.0/4096.0);
//--- allocation
   ArrayResize(state.m_rkc,6);
//--- calculation
   state.m_rkc[0]=37.0/378.0;
   state.m_rkc[1]=0;
   state.m_rkc[2]=250.0/621.0;
   state.m_rkc[3]=125.0/594.0;
   state.m_rkc[4]=0;
   state.m_rkc[5]=512.0/1771.0;
//--- allocation
   ArrayResize(state.m_rkcs,6);
//--- calculation
   state.m_rkcs[0]=2825.0/27648.0;
   state.m_rkcs[1]=0;
   state.m_rkcs[2]=18575.0/48384.0;
   state.m_rkcs[3]=13525.0/55296.0;
   state.m_rkcs[4]=277.0/14336.0;
   state.m_rkcs[5]=1.0/4.0;
   state.m_rkk.Resize(6,n);
//--- Main cycle consists of two iterations:
//--- * outer where we travel from X[i-1] to X[i]
//--- * inner where we travel inside [X[i-1],X[i]]
   state.m_ytbl.Resize(m,n);
   ArrayResize(state.m_escale,n);
   ArrayResize(state.m_yn,n);
   ArrayResize(state.m_yns,n);
//--- change value
   xc=state.m_xg[0];
   for(i_=0; i_<=n-1; i_++)
      state.m_ytbl.Set(0,i_,state.m_yc[i_]);
   for(j=0; j<=n-1; j++)
      state.m_escale[j]=0;
   i=1;
//--- check
   if(i>m-1)
     {
      state.m_repterminationtype=1;
      //--- return result
      return(false);
     }
//--- begin inner iteration
   return(Func_lbl_6(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
  }
//+------------------------------------------------------------------+
//| Auxiliary function for ODESolverIteration. Is a product to get   |
//| rid of the operator unconditional jump goto                      |
//+------------------------------------------------------------------+
void CODESolver::Func_lbl_rcomm(CODESolverState &state,int n,int m,
                                int i,int j,int k,int klimit,
                                bool gridpoint,double xc,double v,
                                double h,double h2,double err,
                                double maxgrowpow)
  {
//--- save
   state.m_rstate.ia.Set(0,n);
   state.m_rstate.ia.Set(1,m);
   state.m_rstate.ia.Set(2,i);
   state.m_rstate.ia.Set(3,j);
   state.m_rstate.ia.Set(4,k);
   state.m_rstate.ia.Set(5,klimit);
   state.m_rstate.ba[0]= gridpoint;
   state.m_rstate.ra.Set(0,xc);
   state.m_rstate.ra.Set(1,v);
   state.m_rstate.ra.Set(2,h);
   state.m_rstate.ra.Set(3,h2);
   state.m_rstate.ra.Set(4,err);
   state.m_rstate.ra.Set(5,maxgrowpow);
  }
//+------------------------------------------------------------------+
//| Auxiliary function for ODESolverIteration. Is a product to get   |
//| rid of the operator unconditional jump goto                      |
//+------------------------------------------------------------------+
bool CODESolver::Func_lbl_6(CODESolverState &state,int &n,int &m,
                            int &i,int &j,int &k,int &klimit,
                            bool &gridpoint,double &xc,double &v,
                            double &h,double &h2,double &err,
                            double &maxgrowpow)
  {
//--- truncate step if needed (beyond right boundary).
//--- determine should we store X or not
   if(xc+h>=state.m_xg[i])
     {
      h=state.m_xg[i]-xc;
      gridpoint=true;
     }
   else
      gridpoint=false;
//--- Update error scale maximums
//--- These maximums are initialized by zeros,
//--- then updated every iterations.
   for(j=0; j<=n-1; j++)
      state.m_escale[j]=MathMax(state.m_escale[j],MathAbs(state.m_yc[j]));
//--- make one step:
//--- 1. calculate all info needed to do step
//--- 2. update errors scale maximums using values/derivatives
//---    obtained during (1)
//--- Take into account that we use scaling of X to reduce task
//--- to the form where x[0] < x[1] < ... < x[n-1]. So X is
//--- replaced by x=xscale*t,and dy/dx=f(y,x) is replaced
//--- by dy/dt=xscale*f(y,xscale*t).
   for(int i_=0; i_<=n-1; i_++)
      state.m_yn[i_]=state.m_yc[i_];
   for(int i_=0; i_<=n-1; i_++)
      state.m_yns[i_]=state.m_yc[i_];
   k=0;
//--- function call, return result
   return(Func_lbl_8(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
  }
//+------------------------------------------------------------------+
//| Auxiliary function for ODESolverIteration. Is a product to get   |
//| rid of the operator unconditional jump goto                      |
//+------------------------------------------------------------------+
bool CODESolver::Func_lbl_8(CODESolverState &state,int &n,int &m,
                            int &i,int &j,int &k,int &klimit,
                            bool &gridpoint,double &xc,double &v,
                            double &h,double &h2,double &err,
                            double &maxgrowpow)
  {
//--- check
   if(k>5)
      return(Func_lbl_10(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
//--- prepare data for the next update of YN/YNS
   state.m_x=state.m_xscale*(xc+state.m_rka[k]*h);
   for(int i_=0; i_<=n-1; i_++)
      state.m_y[i_]=state.m_yc[i_];
//--- calculation
   for(j=0; j<=k-1; j++)
     {
      v=state.m_rkb[k][j];
      for(int i_=0; i_<=n-1; i_++)
         state.m_y[i_]=state.m_y[i_]+v*state.m_rkk[j][i_];
     }
   state.m_needdy=true;
   state.m_rstate.stage=0;
//--- Saving state
   Func_lbl_rcomm(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow);
//--- return result
   return(true);
  }
//+------------------------------------------------------------------+
//| Auxiliary function for ODESolverIteration. Is a product to get   |
//| rid of the operator unconditional jump goto                      |
//+------------------------------------------------------------------+
bool CODESolver::Func_lbl_10(CODESolverState &state,int &n,int &m,
                             int &i,int &j,int &k,int &klimit,
                             bool &gridpoint,double &xc,double &v,
                             double &h,double &h2,double &err,
                             double &maxgrowpow)
  {
//--- estimate error
   err=0;
   for(j=0; j<=n-1; j++)
     {
      //--- check
      if(!state.m_fraceps)
        {
         //--- absolute error is estimated
         err=MathMax(err,MathAbs(state.m_yn[j]-state.m_yns[j]));
        }
      else
        {
         //--- Relative error is estimated
         v=state.m_escale[j];
         //--- check
         if(v==0.0)
            v=1;
         err=MathMax(err,MathAbs(state.m_yn[j]-state.m_yns[j])/v);
        }
     }
//--- calculate new step,restart if necessary
   if(maxgrowpow*err<=state.m_eps)
      h2=m_odesolvermaxgrow*h;
   else
      h2=h*MathPow(state.m_eps/err,0.2);
//--- check
   if(h2<h/m_odesolvermaxshrink)
      h2=h/m_odesolvermaxshrink;
//--- check
   if(err>state.m_eps)
     {
      h=h2;
      //--- begin inner iteration
      return(Func_lbl_6(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
     }
//--- advance position
   xc=xc+h;
   for(int i_=0; i_<=n-1; i_++)
      state.m_yc[i_]=state.m_yn[i_];
//--- update H
   h=h2;
//--- break on grid point
   if(gridpoint)
     {
      //--- save result
      for(int i_=0; i_<=n-1; i_++)
         state.m_ytbl.Set(i,i_,state.m_yc[i_]);
      i=i+1;
      //--- check
      if(i>m-1)
        {
         state.m_repterminationtype=1;
         return(false);
        }
      //--- begin inner iteration
      return(Func_lbl_6(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
     }
//--- begin inner iteration
   return(Func_lbl_6(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
  }
//+------------------------------------------------------------------+
