HDMadapt:=proc(EqODEs,bc1,bc2,pars,xpars,nder,nele,Range,atol,rtol)
    local hnew,N,ics,maxiter,ErrMax,k,sol,ErrList,ErrMin,x0,i,x0old,OKpts,Rmvindx,Wntindx,x0wnt,hwnt,ErrListwnt,hpredict,Nsub,ndernew,
          maxnele,maxnder,maxsize,cnt,rmv,nODEs,Ntot,NN,adaptflag,catchflag,nWntindx,nxpars;
 
    # Copyright: Dr. Venkat R.Subramanian and MAPLE group members at UW.
    # Original code is written by V. R. Subramanian (vsubram@uw.edu).
    # Last modified by Jerry Chen (jerrysyc@uw.edu) on 2017/07/28
 
    optimize;
 
    # Use uniform mesh as a start
    hnew:=[seq((Range[2]-Range[1])/nele,i=1..nele)]:
    N:=nele;
    #ics:=ic0;
    ndernew:=nder;
    nODEs:=nops(EqODEs);
    nxpars:=nops(unknownpars); 
    ics:=[seq(0.5,i=1..(N+1)*nODEs+nxpars)];
  
    # maxiter for adapt
    maxiter:= 20;
    ErrMax:=10;
    maxnele:=200;
    maxnder:=9;
    maxsize:=500; # maximum size of system
 
   
    # loop till we find the converged answer  
    for k to maxiter while ErrMax > atol or ndernew < maxnder+1 do
  
      try
 
       sol:=HDM(EqODEs,bc1,bc2,pars,xpars,ndernew,N,Range,ics,hnew,atol,rtol);
 
       # store the norm error as a list
       ErrList:=op(sol[2]);
  
       # find the min and max for the judging criteria
       ErrMax:=max(ErrList);
       ErrMin:=min(ErrList);
 
       catchflag:= 0;
 
     catch:
       print("HDM fail, increase points to...", 2*N);
       N:=2*N;
       hnew:=[seq((Range[2]-Range[1])/N,i=1..N)]:
       ics:=[seq(0.5,i=1..(N+1)*nODEs+nxpars)];
       if ndernew <= maxnder and N > maxnele then ndernew:= ndernew+1; end if;
       catchflag:= 1;
     end try:
 
      # confrim to find the converged answer
      if ErrMax < atol then break; end;
   
      # if the situation is really bad (error is really big) or cannot find the answer => doulbe the nodes
      if ErrMin > 0.1 and catchflag = 0 then 
        print("ErrMin > 0.1, double the mesh to...",2*N);
        hnew:=[seq(seq(hnew[i]/2,j=1..2),i=1..N)];
        NN:=[seq(op([1,seq(2,i=1..N)]),j=1..nODEs)];
        ics:=[seq(seq(seq(map(rhs,sol[1])[(k-1)*(N+1)+i],j=1..NN[i]),i=1..N+1),k=1..nODEs),
              seq(rhs(sol[1][nODEs*(N+1)+i]),i=1..nxpars)];
        N:=2*N;   
  
      # minimize the error by insertion local points
      elif catchflag = 0 then
  
        # initial mesh points before add and removal
        x0:=Vector(N+1,datatype=float[8]);
        x0[1]:=Range[1]:
        for i from 2 to N+1 do x0[i]:=x0[i-1]+hnew[i-1]:od: 
        x0old:=[seq(x0[i],i=1..N+1)]; #print("initial mesh points before add and removal",x0old);
        
        # Find the index of the norm ErrList which is less then (atol/10)          
        OKpts:=[seq(`if`(ErrList[i]<atol/10,i,NULL),i=1..N)];
  
        # based on OKpts, find the which index should be removed
        # Removal rules: in a grouping of 3, remove 2 (4) point from 1,2,3,(4,5)
        cnt:=1;
        for i from 1 to nops(OKpts)-2 do
          if OKpts[i+2]=OKpts[i]+2 then
            rmv[cnt]:=i+1;
            i:=i+1;
            cnt:=cnt+1;
          end if;
        end do;
        Rmvindx:=[seq(rmv[i],i=1..cnt-1)]; 
        Wntindx:=remove(has,[seq(i,i=1..N)],Rmvindx);
        nWntindx:= nops(Wntindx);
     
        # New nodes and h after removal
        x0wnt:=[x0[1],seq(x0[i+1],i=Wntindx)];
        hwnt:=[seq(x0wnt[i+1]-x0wnt[i],i=1..nWntindx)];
        ErrListwnt:=[seq(ErrList[i],i=Wntindx)];
  
        # Find the new h based on the 2n-1 accuracy
        ## Nsub = number of insertion points for subintervals
        ## hpredict = predict spacing
        for i from 1 to nWntindx do hpredict[i]:=hwnt[i]*(atol/ErrListwnt[i])^(1/(2*nder-1));od;
        for i from 1 to nWntindx do Nsub[i]:=trunc(hwnt[i]/hpredict[i])+1; od;
        Ntot:=add(Nsub[i],i=1..nWntindx);    
  
        # update the info for the new mesh
        if (Ntot > maxnele or Ntot*nODEs > maxsize) or adaptflag = 1 then 
          print("case1 - double the nodes to...",2*N);
          hnew:=[seq(seq(hnew[i]/2,j=1..2),i=1..N)];
          NN:=[seq(op([1,seq(2,i=1..N)]),j=1..nODEs)];
          if adaptflag = 1 then
            ics:=[seq(seq(seq(map(rhs,sol[1])[(k-1)*(N+1)+i],j=1..NN[i]),i=1..N+1),k=1..nODEs),
                  seq(rhs(sol[1][nODEs*(N+1)+i]),i=1..nxpars)];
          else
            ics:=[seq(0.5,i=1..(2*N+1)*nODEs+nxpars)];
          end if:
          if ndernew <= maxnder and N > maxnele then ndernew:= ndernew+1; end if:
          N:=2*N; 
          adaptflag:= 0;
        else 
          print("case2 - insert/remove points based on 2n-1 error to...",Ntot);
          hnew:=[seq(seq(hwnt[i]/Nsub[i],j=1..Nsub[i]),i=1..nWntindx)];
          NN:=[seq(op([1,seq(Nsub[i],i=1..nWntindx)]),j=1..nODEs)]; 
          ics:=[seq(seq(seq(map(rhs,sol[1])[(k-1)*(nWntindx+1)+i],j=1..NN[i]),i=1..nWntindx+1),k=1..nODEs),
                seq(rhs(sol[1][nODEs*(nWntindx+1)+i]),i=1..nxpars)];
          N:=Ntot;   
          adaptflag:= 1;   
        end if;
        
      end if; # End with decision if ErrMin > 0.1  
  
    end do; # close the loop for maxiter
  
    # Return 
    [sol[1],([seq(evalf(ErrList[i]),i=1..N)]),hnew,ndernew,ErrMax];
  
 end proc:

 HDM:=proc(eqodes,bc1,bc2,pars,xpars,nder,N,rng,ic,hh,atol,rtol)
    local node,odevars,vars,f,i,F,x0,pade_alpha,pade_beta,alpha,beta,alpha1,beta1,YY,Eqq,Eqs,varsNR,iter,guess,sol,Err,nxpars;
  
    # Copyright: Dr. Venkat R.Subramanian and MAPLE group members at UW.
    # Original code is written by V. R. Subramanian (vsubram@uw.edu).
    # Last modified by Jerry Chen (jerrysyc@uw.edu) on 2017/08/03
   
    optimize;
  
    # Define the variables 
    node:=nops(eqodes);
    odevars:=select(type,map(op,map(lhs,eqodes)),'function');
    vars:=[seq(odevars[i]=Y||i[j],i=1..node)]; # i=vars, j=nodes
 
    # Find the derivatives 
    f[1]:=Vector(node,[seq(rhs(eqodes[i]),i=1..node)]);
 
    for i from 2 to nder do 
      f[i]:=simplify(subs(eqodes,Vector(node,[seq(diff(f[i-1][j],x),j=1..node)])));
    end do:
 
    for i from 1 to nder do 
      F[i]:=unapply(subs(op(pars),vars,f[i]),j,x);
    od:
 
    # Initialize the mesh
    x0:=Vector(N+1,datatype=float[8]):  
    x0[1]:=rng[1]:
    for i from 1 to N do x0[i+1]:=x0[i]+hh[i]: od:
  
    # Pade coefficients
    pade_alpha:=unapply((m+n-ii)!*(m)!/((m+n)!*ii!*(m-ii)!),m,n,ii):
    pade_beta:=unapply((m+n-ii)!*(n)!/((m+n)!*ii!*(n-ii)!),m,n,ii):
 
    for i from 1 to nder do
      alpha[i]:= pade_alpha(nder,nder,i);
      beta[i]:= alpha[i];
      beta1[i]:= pade_beta(nder-1,nder,i); 
    end do:
 
    for i from 1 to nder-1 do
      alpha1[i]:=pade_alpha(nder-1,nder,i); 
    end do:
    alpha1[nder]:=0;
 
    # HDM formula
    YY:=unapply(Vector(node,[seq(Y||i[j],i=1..node)]),j);
    for i from 1 to N do 
      Eqq[i]:=eval(YY(i-1)-YY(i)+add(hh[i]^j*(alpha[j]*F[j](i-1,x0[i])+(-1)^(j+1)*beta[j]*F[j](i,x0[i+1])),j=1..nder));
    od;
    
    Eqs:=evalf([seq(seq(Eqq[i][j],i=1..N),j=1..node),op(subs(op(pars),vars,j=0,x=rng[1],bc1)),op(subs(op(pars),vars,j=N,x=rng[2],bc2))]);
 
    # vars for NR
    varsNR:=[seq(seq(Y||j[i],i=0..N),j=1..node),op(xpars)]:
  
    # NR
    iter:=15;
    nxpars:=nops(xpars);
    guess:=[seq(varsNR[i]=ic[i],i=1..node*(N+1)+nxpars)];
    sol:=SNewton(Eqs,guess,atol,rtol,iter,N,nxpars): 
  
    # Error (2n HDM formula - (2n-1) HDM formula) - Note: not consider error at the first point, x0
    for i from 1 to N do 
      Err[i]:=evalf(LinearAlgebra:-Norm(subs(sol,add(hh[i]^j*((alpha[j]-alpha1[j])*F[j](i-1,x0[i])
              +(-1)^(j+1)*(beta[j]-beta1[j])*F[j](i,x0[i+1])),j=1..nder)))/sqrt(node));
    end do;
  
    # Return 
    [sol,([seq(evalf(Err[i]),i=1..N)])];
  
 end proc:

 SNewton := proc (Eqs, ics, atol, rtol, iter, nele, nxpars) 
     local Jac,N,scale,vars,Eqs1,yold,jac,i,LL,j,F,k,errL,dy,ynew,rhsics,N1,bcindex; 
     
     # Copyright: Dr. Venkat R.Subramanian and MAPLE group members at UW.
     # Original code is written by V. R. Subramanian (vsubram@uw.edu).
     # Last modified by Jerry Chen (jerrysyc@uw.edu) on 2017/08/03
  
     optimize;
     
     # Initialization
     N:=nops(Eqs);
     errL := 10; 
     
     # Input check
     if nops(Eqs) <> nops(ics) then ERROR("Check the number of equations and variables"); end;
     if nops(remove(has,indets(Eqs,name),map(lhs,ics))) <> 0 then ERROR("Variables' name don't match"); end;
  
     # Scale based on the ics to force all y calculate from 0 to 1.
     rhsics:=map(rhs,ics);
     vars:=[seq(lhs(ics[i])=y[i],i=1..N)];
     yold :=evalf(Vector([seq(rhsics[i],i=1..N)]));
 
     # Create equation set
     Eqs1:=subs(vars,Eqs);
     F := unapply(`<,>`(seq(Eqs1[i],i = 1 .. N)),y); 
   
     # Create Jacobian
     jac := Matrix(1 .. N,1 .. N,storage = sparse);
     N1:=(N-nxpars)/(nele+1); 
     for k from 1 to N1 do
       for i from 1 to nele do
         for j from 1 to N1 do
           jac[i+(k-1)*nele,i+(j-1)*(nele+1)]:=diff(Eqs1[i+(k-1)*(nele)],y[i+(j-1)*(nele+1)]); 
           jac[i+(k-1)*nele,i+1+(j-1)*(nele+1)]:=diff(Eqs1[i+(k-1)*(nele)],y[i+1+(j-1)*(nele+1)]);
         od:
       od:
     od:
     bcindex:=[seq(1+(i-1)*(nele+1),i=1..N1),seq(nele+1+(i-1)*(nele+1),i=1..N1)];
     for i from N1*nele+1 to N-nxpars do 
       for j from 1 to nops(bcindex) do
         jac[i,bcindex[j]]:=diff(Eqs1[i],y[bcindex[j]]);
       od:
     od:    
     for j from N-nxpars+1 to N do
       for i from 1 to N do
         jac[i,j]:=diff(Eqs1[i],y[j]):
         jac[j,i]:=diff(Eqs1[j],y[i]):
       od:
     od:
     Jac:=unapply(jac,y);
  
     # Iterate until the y converges less than the tolerance
     for k to iter while errL > 1 do
  
       # Avoid version caused "SparseDirect" failed happens in Maple classic worksheet 
       dy := LinearAlgebra:-LinearSolve(-Jac(yold),F(yold),method=SparseDirect); 
 
       if add(`if`(type((dy[i]),numeric),0,1),i=1..N)>0 then ERROR("Cannot converge"); end if;
       ynew := yold + dy; 
       errL := LinearAlgebra:-Norm(dy,2)/LinearAlgebra:-Norm(Vector([seq(atol+rtol*ynew[i],i=1..N)]),2);
 
       # update yold to next iteration 
       yold := ynew; 
  
     end do; 
    
     if iter < k then ERROR("Exceed max NR iteration: %1",iter) end if; 
  
     # Return
     [seq(lhs(ics[i])=ynew[i],i=1..N)]; # with name
 
     
 end proc:

 HDMpars:=proc(eqodes,bc1,bc2,pars,xpars,nder,nele,rng,hh,atol,rtol)
     local node,odevars,vars,f,i,F,x0,pade_alpha,pade_beta,alpha,beta,alpha1,beta1,YY,Eqq,Eqs,varsNR,guess,sol,Err,nxpars,
           N,errL,Eqs1,jac,N1,k,j,bcindex,Jac,LErr;
   
     # Copyright: Dr. Venkat R.Subramanian and MAPLE group members at UW.
     # Original code is written by V. R. Subramanian (vsubram@uw.edu).
     # Last modified by Jerry Chen (jerrysyc@uw.edu) on 2017/08/11
    
     optimize;
   
     # Define the variables 
     node:=nops(eqodes);
     odevars:=select(type,map(op,map(lhs,EqODEs)),'function');
     vars:=[seq(odevars[i]=Y||i[j],i=1..node)]; # i=vars, j=nodes
  
     # Find the derivatives 
     f[1]:=Vector(node,[seq(rhs(eqodes[i]),i=1..node)]);
  
     for i from 2 to nder do 
       f[i]:=simplify(subs(eqodes,Vector(node,[seq(diff(f[i-1][j],x),j=1..node)])));
     end do:
  
     for i from 1 to nder do 
       F[i]:=unapply(subs(vars,f[i]),j,x);
     od:
  
     # Initialize the mesh
     x0:=Vector(nele+1,datatype=float[8]):  
     x0[1]:=rng[1]:
     for i from 1 to nele do x0[i+1]:=x0[i]+hh[i]: od:
   
     # Pade coefficients
     pade_alpha:=unapply((m+n-ii)!*(m)!/((m+n)!*ii!*(m-ii)!),m,n,ii):
     pade_beta:=unapply((m+n-ii)!*(n)!/((m+n)!*ii!*(n-ii)!),m,n,ii):
  
     for i from 1 to nder do
       alpha[i]:= pade_alpha(nder,nder,i);
       beta[i]:= alpha[i];
       beta1[i]:= pade_beta(nder-1,nder,i); 
     end do:
  
     for i from 1 to nder-1 do
       alpha1[i]:=pade_alpha(nder-1,nder,i); 
     end do:
     alpha1[nder]:=0;
  
     # HDM formula
     YY:=unapply(Vector(node,[seq(Y||i[j],i=1..node)]),j);
     for i from 1 to nele do 
       Eqq[i]:=eval(YY(i-1)-YY(i)+add(hh[i]^j*(alpha[j]*F[j](i-1,x0[i])+(-1)^(j+1)*beta[j]*F[j](i,x0[i+1])),j=1..nder));
     od;
     
     Eqs:=evalf([seq(seq(Eqq[i][j],i=1..nele),j=1..node),op(subs(vars,j=0,x=rng[1],bc1)),op(subs(vars,j=nele,x=rng[2],bc2))]);
 
     # vars for NR
     varsNR:=[seq(seq(Y||j[i],i=0..nele),j=1..node),op(xpars)]:
  
     # NR
     nxpars:=nops(xpars);
     N:=nops(Eqs);
     vars:=[seq(varsNR[i]=y[i],i=1..N)];
  
     # Local Error
     LErr:=Vector([seq(evalf(LinearAlgebra:-Norm(subs(vars,add(hh[i]^j*((alpha[j]-alpha1[j])*F[j](i-1,x0[i])
               +(-1)^(j+1)*(beta[j]-beta1[j])*F[j](i,x0[i+1])),j=1..nder)))/sqrt(node)),i=1..nele)]);
 
     # Create equation set
     Eqs1:=subs(vars,Eqs);
      
     # Create Jacobian
     jac := Matrix(1 .. N,1 .. N,storage = sparse);
     N1:=(N-nxpars)/(nele+1); 
     for k from 1 to N1 do
       for i from 1 to nele do
         for j from 1 to N1 do
           jac[i+(k-1)*nele,i+(j-1)*(nele+1)]:=diff(Eqs1[i+(k-1)*(nele)],y[i+(j-1)*(nele+1)]); 
           jac[i+(k-1)*nele,i+1+(j-1)*(nele+1)]:=diff(Eqs1[i+(k-1)*(nele)],y[i+1+(j-1)*(nele+1)]);
         od:
       od:
     od:
     bcindex:=[seq(1+(i-1)*(nele+1),i=1..N1),seq(nele+1+(i-1)*(nele+1),i=1..N1)];
     for i from N1*nele+1 to N-nxpars do 
       for j from 1 to nops(bcindex) do
         jac[i,bcindex[j]]:=diff(Eqs1[i],y[bcindex[j]]);
       od:
     od:    
     for j from N-nxpars+1 to N do
       for i from 1 to N do
         jac[i,j]:=diff(Eqs1[i],y[j]):
         jac[j,i]:=diff(Eqs1[j],y[i]):
       od:
     od:
        
     [Eqs1,jac,LErr,node];
   
 end proc:


 HDMpars2 := proc (Gen,ic0,pars) 
      local Jac,N,scale,vars,Eqs1,yold,jac,i,LL,j,F,k,errL,dy,ynew,rhsics,N1,bcindex,iter,Lerr,node; 
      
      # Copyright: Dr. Venkat R.Subramanian and MAPLE group members at UW.
      # Original code is written by V. R. Subramanian (vsubram@uw.edu).
      # Last modified by Jerry Chen (jerrysyc@uw.edu) on 2017/08/11
   
      optimize;
 
      # Initialization
      N:=nops(Gen[1]);
      errL := 10; 
      iter:=15;
      node:=Gen[4];
 
      F := unapply(`<,>`(seq(subs(op(pars),Gen[1][i]),i = 1 .. N)),y); 
      Jac:=unapply(subs(op(pars),Gen[2]),y);
 
      # Scale based on the ics to force all y caculate from 0 to 1.
      rhsics:=ic0;   #map(rhs,ics);
      yold:=evalf(Vector([seq(rhsics[i],i=1..N)]));
      
      # Iterate until the y converges less then the tolerence
      for k to iter while errL > 1 do
   
        # Avoid version caused "SparseDirect" failed happens in Maple classic worksheet 
        dy := LinearAlgebra:-LinearSolve(-Jac(yold),F(yold),method=SparseDirect); 
  
        if add(`if`(type((dy[i]),numeric),0,1),i=1..N)>0 then ERROR("Cannot converge"); end if;
        ynew := yold + dy; 
        errL := LinearAlgebra:-Norm(dy,2)/LinearAlgebra:-Norm(Vector([seq(atol+rtol*ynew[i],i=1..N)]),2);
  
        # update yold to next iteration 
        yold := ynew; 
   
      end do; 
     
      if iter < k then ERROR("Exceed max NR iteration: %1",iter) end if; 
   
      Lerr:=subs(seq(y[i]=ynew[i],i=1..N),op(pars),Gen[3]);
        
      # Return
      [ynew,Lerr]; # without name
    
 end proc: