#Code written by Venkat R. Subramanian and MAPLE group members, UT Austin, last updated on 11/30/2022. Released under MIT license
LUDecomp:=define_external("hw_SpUMFPACK_MatFactor",':-MAPLE', ':-LIB' = "linalg"):
LUSolve:=define_external("hw_SpUMFPACK_MatMatSolve",':-MAPLE', ':-LIB' = "linalg"):
FreeMat:=define_external("hw_SpUMFPACK_FreeNumeric",':-MAPLE', ':-LIB' = "linalg"):
  

IMPDAE := proc(Eqode, Eqae1, ICs,solveroptions)

local node, nae,  j, eqs,Vars1, varsold, n1,count,
varsdir, i, ss, Equation, k, F, Vars, j11, L, LL, LL2,
Jac, ic0, y0, yold, err, dt, singlestep, tol, Yf,t11,Nt,iter1,icb,
dt1,Anz,NumericA,J1,ymid,yh,yh2,yR,tt,ee,atol1,rtol1,errL,growth,nextH,hold,kmod1,i1,tcont,YY,
tf,atol,hinit,hmax,Ntot,iter,kmod,kjac,fstep,varsdir1,solf;
optimize;
with(LinearAlgebra):

tf:=rhs(solveroptions[1]):atol:=rhs(solveroptions[2]):hinit:=rhs(solveroptions[3]):hmax:=rhs(solveroptions[4]):
Ntot:=rhs(solveroptions[5]):iter:=rhs(solveroptions[6]):


YY:=Array(1..Ntot+1):
    
    tcont:=1;
    node := nops(Eqode);
    nae := nops(Eqae1);
    Nt:=node+nae:
    Vars := [seq(lhs(ICs[i]), i = 1 .. Nt)];
   kjac:=0:
   fstep:=0: 
yh:=Vector(Nt,datatype=float[8]):
yh2:=Vector(Nt,datatype=float[8]):
ymid:=Vector(Nt,datatype=float[8]):
eqs:=Array(1..node+nae):
for i from 1 to node do eqs[i]:=rhs(Eqode[i]):od:
for i from 1 to nae do eqs[i+node]:=lhs(Eqae1[i]) - rhs(Eqae1[i]):od:
Equation:=Array(1..node+nae):
j11 := Matrix(1 .. Nt, 1 .. Nt, storage = sparse);
J1:=Matrix(1 .. Nt, 1 .. Nt, storage = sparse);
for i from 1 to Nt do 
L:=indets(eqs[i])minus{t};
LL := [seq(ListTools:-Search(L[j], Vars),
             j = 1 .. nops(L))];
 LL2 := ListTools:-MakeUnique(LL);
if LL2[nops(LL2)] = 0 then
             LL := [seq(LL2[i1], i1 = 1 .. nops(LL2) - 1)]
         end if;
         if LL2[1]= 0 then LL:=[seq(LL2[i1],i1=2..nops(LL2))]: end:#print(LL);
varsdir:=[seq(Vars[LL[i]]=uu[LL[i]]/2+uu_old[LL[i]],i=1..nops(LL))];
varsdir1:=[seq(Vars[LL[i]]=uu[LL[i]]+uu_old[LL[i]],i=1..nops(LL))];
varsold:=[seq(Vars[LL[i]]=uu_old[LL[i]],i=1..nops(LL))];
if i<=node then Equation[i]:=uu[i]-deltA*subs(varsdir,eqs[i])-deltA*0*subs(varsold,eqs[i]):

else Equation[i]:=subs(varsdir1,eqs[i]):

end:
od:
Vars1 := [seq(uu[i], i = 1 .. Nt)];
    
    for i to Nt do
        
        L := indets(Equation[i]);
        LL := [seq(ListTools:-Search(L[j], Vars1),
            j = 1 .. nops(L))];
        LL2 := ListTools:-MakeUnique(LL);
        if LL2[nops(LL2)] = 0 then
            LL := [seq(LL2[i], i = 1 .. nops(LL2) - 1)]
        end if;
        if LL2[1]= 0 then LL:=[seq(LL2[i],i=2..nops(LL2))]: end:

        for j to nops(LL) do
            #if LL[j]>=i then j11[i, LL[j]] := diff(Eqss[i], uu[LL[j]]);end:
            #j11[i, LL[j]] := diff(Eqss[i], uu[LL[j]]);
            j11[i, LL[j]] := diff(Equation[i], uu[LL[j]]);
        end do;
    end do;

  
F := unapply(<seq(evalf(Equation[i]), i = 1 .. node + nae)>, uu,
        uu_old, deltA);#print(F);

    Jac := unapply(j11, uu, uu_old, deltA);#print(Jac);
    ic0 := Vector(node + nae,
        [seq(rhs(ICs[i]), i = 1 .. node + nae)],
        datatype = float[8]);
    y0 := Vector(Nt, datatype = float[8]);
    yold := y0;
    err := 1;
    
    singlestep := proc(ic0,y00, h, node, nae,iter,Anz,NumericA)
global ncons;
        local k, db, ynew, err, yold, i, dt, ic1, ydb,Nt,J1,F1;
            Nt := node + nae;
            yold:=copy(y00):
            db:=copy(y00):
            err := 1;
            
            for k to iter while 1e-14 < err do
            
            F1:=-convert(F(yold, ic0, h),Matrix,storage=sparse):
            db:=convert(LUSolve(Nt,1,Anz,NumericA,F1),Vector):
                ynew := yold + db;
                err := evalhf(eval(LinearAlgebra:-Norm(db, Nt)/LinearAlgebra:-Norm(ynew, Nt)));
                yold := ynew;
            end do;
            
            ic1 := ic0 + ynew;
            end proc;
    tol := atol;atol1:=atol:rtol1:=10*atol1:
            count:=0;
            
            
  		yold := Vector(Nt, 0, datatype = float[8]):
                
            try
            count:=count+1:
            tt[count]:=0.0:t11:=tt[count]: 
            J1:= convert(Jac(yold, ic0,0.0),Matrix,datatype=float[8]):
            (Anz,NumericA):=LUDecomp(Nt,Nt,J1):
           
            YY[count] := singlestep(ic0,yold, 0.0, node, nae,11,Anz,NumericA);#print(YY[count],t11):
            catch:
            tcont:=0:
            end try:
		
            dt:=hinit:
            kjac:=0:
            kmod1:=kmod:
            while t11<tf and count<=Ntot and tcont=1 do
            try
            count:=count+1: 
            tt[count]:=tt[count-1]+dt:
            t11:=tt[count]:
            if kjac = 0 then 
            J1:= convert(Jac(yold, YY[count-1], dt),Matrix,datatype=float[8]):
            (Anz,NumericA):=LUDecomp(Nt,Nt,J1):end:
            #yh:=singlestep(YY[count-1],(YY[count-1]-YY[count-2])*dt/hold, dt, node, nae,iter,Anz,NumericA):#print(yh);
            if count=2 then yh:=singlestep(YY[count-1],yold, dt, node, nae,iter,Anz,NumericA):
            else
            yh:=singlestep(YY[count-1],(YY[count-1]-YY[count-2])*dt/hold, dt, node, nae,iter,Anz,NumericA):
            end:
            ymid:=singlestep(YY[count-1], (yh-YY[count-1])/2,dt/2.0, node, nae,iter,Anz,NumericA):
             yh2:=singlestep(ymid,yh-ymid, dt/2.0, node, nae,iter,Anz,NumericA):
            YY[count]:=4.0/3*yh2-yh/3:#print(count,YY[count]): 

            
            
            ee := Vector(node + nae, [seq(abs(yh[i] - yh2[i])
                /(atol1 + rtol1*abs(YY[count][i])),
                i = 1 .. node + nae)]);
                errL := max(atol1,
                1/3*LinearAlgebra:-Norm(ee, 2)/(node + nae));
            if errL<1.0 then 
            growth := min(3.0,max(0.3,0.9*(1/errL)^(1/(2+1))));
            hold:=dt:
            dt:=min(tf-t11,hmax,growth*dt):#print(dt):
            
            
            if errL>0.1 then kjac:=0: else kjac:=1:end:
            else dt:=dt/4: fstep:=fstep+1:#print("step failure at",t11):
            kjac:=0:
            count:=count-1:
            if count=1 then t11:=0: else t11:=tt[count-1]:end:
            if dt<hinit/1000 then tcont:=0:print("cannot integrate past",t11)end:
            end:
            catch:
            
            tcont:=0:
            end try:
            end:
            
if tt[count]=0 then print("Simulation was stopped due to initialization failure"):tcont:=1:end:
if tt[count]<tf and tcont=0 then print("Simulation was stopped before the final time due to a state becoming complex"):end:           
if tt[count]<tf and tcont=1 then print("Simulation was stopped before the final time, increase total time steps or change other solver options"):end: 
if fstep >0 and tcont=1 then print("Integration completed, but number of failed steps=",fstep):end:

            
#[[seq([tt[i1],YY[i1]],i1=1..count)],count];
solf:=Vector(count):
for i from 1 to count do solf[i]:=[t=tt[i],seq(Vars[j]=YY[i][j],j=1..Nt)]:od:
solf;

end proc;


CNDAE := proc(Eqode, Eqae1, ICs,solveroptions)
global ncons;
local node, nae,  j, eqs,Vars1, varsold, n1,count,
varsdir, i, ss, Equation, k, F, Vars, j11, L, LL, LL2,
Jac, ic0, y0, yold, err, dt, singlestep, tol, Yf,t11,Nt,iter1,icb,
dt1,Anz,NumericA,J1,ymid,yh,yh2,yR,tt,ee,atol1,rtol1,errL,growth,nextH,hold,kmod1,i1,tcont,YY,
tf,atol,hinit,hmax,Ntot,iter,kmod,kjac,fstep,varsdir1,solf;
optimize;
with(LinearAlgebra):

tf:=rhs(solveroptions[1]):atol:=rhs(solveroptions[2]):hinit:=rhs(solveroptions[3]):hmax:=rhs(solveroptions[4]):
Ntot:=rhs(solveroptions[5]):iter:=rhs(solveroptions[6]):


YY:=Array(1..Ntot+1):
    tcont:=1;
    node := nops(Eqode);
    nae := nops(Eqae1);
    Nt:=node+nae:
    Vars := [seq(lhs(ICs[i]), i = 1 .. Nt)];
   kjac:=0:
   fstep:=0: 
yh:=Vector(Nt,datatype=float[8]):
yh2:=Vector(Nt,datatype=float[8]):
ymid:=Vector(Nt,datatype=float[8]):
eqs:=Array(1..node+nae):
for i from 1 to node do eqs[i]:=rhs(Eqode[i]):od:
for i from 1 to nae do eqs[i+node]:=lhs(Eqae1[i]) - rhs(Eqae1[i]):od:
Equation:=Array(1..node+nae):
j11 := Matrix(1 .. Nt, 1 .. Nt, storage = sparse);
J1:=Matrix(1 .. Nt, 1 .. Nt, storage = sparse);
for i from 1 to Nt do 
L:=indets(eqs[i])minus{t};
LL := [seq(ListTools:-Search(L[j], Vars),
             j = 1 .. nops(L))];
 LL2 := ListTools:-MakeUnique(LL);
if LL2[nops(LL2)] = 0 then
             LL := [seq(LL2[i1], i1 = 1 .. nops(LL2) - 1)]
         end if;
         if LL2[1]= 0 then LL:=[seq(LL2[i1],i1=2..nops(LL2))]: end:#print(LL);
varsdir:=[seq(Vars[LL[i]]=uu[LL[i]]+uu_old[LL[i]],i=1..nops(LL))];
varsold:=[seq(Vars[LL[i]]=uu_old[LL[i]],i=1..nops(LL))];
if i<=node then Equation[i]:=uu[i]-deltA/2*subs(varsdir,eqs[i])-deltA/2*subs(varsold,eqs[i]):
else Equation[i]:=subs(varsdir,eqs[i]):


end:
od:
Vars1 := [seq(uu[i], i = 1 .. Nt)];
    
    for i to Nt do
        L := indets(Equation[i]);
        LL := [seq(ListTools:-Search(L[j], Vars1),
            j = 1 .. nops(L))];
        LL2 := ListTools:-MakeUnique(LL);
        if LL2[nops(LL2)] = 0 then
            LL := [seq(LL2[i], i = 1 .. nops(LL2) - 1)]
        end if;
        if LL2[1]= 0 then LL:=[seq(LL2[i],i=2..nops(LL2))]: end:

        for j to nops(LL) do
            #if LL[j]>=i then j11[i, LL[j]] := diff(Eqss[i], uu[LL[j]]);end:
            #j11[i, LL[j]] := diff(Eqss[i], uu[LL[j]]);
            j11[i, LL[j]] := diff(Equation[i], uu[LL[j]]);
        end do;
    end do;

  
F := unapply(<seq(evalf(Equation[i]), i = 1 .. node + nae)>, uu,
        uu_old, deltA);#print(F);

    Jac := unapply(j11, uu, uu_old, deltA);#print(Jac);
    ic0 := Vector(node + nae,
        [seq(rhs(ICs[i]), i = 1 .. node + nae)],
        datatype = float[8]);
    y0 := Vector(Nt, datatype = float[8]);
    yold := y0;
    err := 1;
    singlestep := proc(ic0,y00, h, node, nae,iter,Anz,NumericA)
global ncons;
        local k, db, ynew, err, yold, i, dt, ic1, ydb,Nt,J1,F1;
            Nt := node + nae;
            yold:=copy(y00):
            db:=copy(y00):
            err := 1;
            
            for k to iter while 1e-14 < err do
            
            F1:=-convert(F(yold, ic0, h),Matrix,storage=sparse):
            db:=convert(LUSolve(Nt,1,Anz,NumericA,F1),Vector):
                ynew := yold + db;
                err := evalhf(eval(LinearAlgebra:-Norm(db, Nt)/LinearAlgebra:-Norm(ynew, Nt)));
                yold := ynew;
            end do;
            
            ic1 := ic0 + ynew;
            end proc;
    tol := atol;atol1:=atol:rtol1:=10*atol1:
            count:=0;
            
            
  		yold := Vector(Nt, 0, datatype = float[8]):
               
            try
            count:=count+1:
            tt[count]:=0.0:t11:=tt[count]: 
            J1:= convert(Jac(yold, ic0,0.0),Matrix,datatype=float[8]):
            (Anz,NumericA):=LUDecomp(Nt,Nt,J1):
           
            YY[count] := singlestep(ic0,yold, 0.0, node, nae,11,Anz,NumericA);#print(YY[count],t11):
            catch:
            tcont:=0:
            end try:
		
            dt:=hinit:
            kjac:=0:
            
            while t11<tf and count<=Ntot and tcont=1 do
            try
            count:=count+1: 
            tt[count]:=tt[count-1]+dt:
            t11:=tt[count]:
            
            if kjac = 0 then 
            J1:= convert(Jac(yold, YY[count-1], dt),Matrix,datatype=float[8]):
            (Anz,NumericA):=LUDecomp(Nt,Nt,J1):end:
            #yh:=singlestep(YY[count-1],yold, dt, node, nae,iter,Anz,NumericA):#print(yh);
            if count=2 then yh:=singlestep(YY[count-1],yold, dt, node, nae,iter,Anz,NumericA):
            else
            yh:=singlestep(YY[count-1],(YY[count-1]-YY[count-2])*dt/hold, dt, node, nae,iter,Anz,NumericA):
            end:
            ymid:=singlestep(YY[count-1], (yh-YY[count-1])/2,dt/2.0, node, nae,iter,Anz,NumericA):
             yh2:=singlestep(ymid,yh-ymid, dt/2.0, node, nae,iter,Anz,NumericA):
            YY[count]:=4.0/3*yh2-yh/3:#print(count,YY[count]): 

            
            
            ee := Vector(node + nae, [seq(abs(yh[i] - yh2[i])
                /(atol1 + rtol1*abs(YY[count][i])),
                i = 1 .. node + nae)]);
                errL := max(atol1,
                1/3*LinearAlgebra:-Norm(ee, 2)/(node + nae));
            if errL<1.0 then 
            growth := min(3.0,max(0.3,0.9*(1/errL)^(1/(2+1))));
            hold:=dt:
            dt:=min(tf-t11,hmax,growth*dt):#print(dt):
            
            
            if errL>0.1 then kjac:=0: else kjac:=1:end:
            else dt:=dt/4: fstep:=fstep+1:#print("step failure at",t11):
            kjac:=0:
            count:=count-1:
             if count=1 then t11:=0: else t11:=tt[count-1]:end:
            if dt<hinit/1000 then tcont:=0:print("cannot integrate past",t11)end:
            end:
            catch:
            
            tcont:=0:
            end try:
            end:
            
if tt[count]=0 then print("Simulation was stopped due to initialization failure"):tcont:=1:end:
if tt[count]<tf and tcont=0 then print("Simulation was stopped before the final time due to a state becoming complex"):end:           
if tt[count]<tf and tcont=1 then print("Simulation was stopped before the final time, increase total time steps or change other solver options"):end: 
if fstep >0 and tcont=1 then print("Integration completed, but number of failed steps=",fstep):end:
            
#[[seq([tt[i1],YY[i1]],i1=1..count)],count];
solf:=Vector(count):
for i from 1 to count do solf[i]:=[t=tt[i],seq(Vars[j]=YY[i][j],j=1..Nt)]:od:
solf;

end proc;

EBDAE := proc(Eqode, Eqae1, ICs,solveroptions)
global ncons;
local node, nae,  j, eqs,Vars1, varsold, n1,count,
varsdir, i, ss, Equation, k, F, Vars, j11, L, LL, LL2,
Jac, ic0, y0, yold, err, dt, singlestep, tol, Yf,t11,Nt,iter1,icb,
dt1,Anz,NumericA,J1,ymid,yh,yh2,yR,tt,ee,atol1,rtol1,errL,growth,nextH,hold,kmod1,i1,tcont,YY,
tf,atol,hinit,hmax,Ntot,iter,kmod,kjac,fstep,varsdir1,solf;
optimize;
with(LinearAlgebra):

tf:=rhs(solveroptions[1]):atol:=rhs(solveroptions[2]):hinit:=rhs(solveroptions[3]):hmax:=rhs(solveroptions[4]):
Ntot:=rhs(solveroptions[5]):iter:=rhs(solveroptions[6]):


YY:=Array(1..Ntot+1):
    tcont:=1;
    node := nops(Eqode);
    nae := nops(Eqae1);
    Nt:=node+nae:
    Vars := [seq(lhs(ICs[i]), i = 1 .. Nt)];
   kjac:=0:
   fstep:=0: 
yh:=Vector(Nt,datatype=float[8]):
yh2:=Vector(Nt,datatype=float[8]):
ymid:=Vector(Nt,datatype=float[8]):
eqs:=Array(1..node+nae):
for i from 1 to node do eqs[i]:=rhs(Eqode[i]):od:
for i from 1 to nae do eqs[i+node]:=lhs(Eqae1[i]) - rhs(Eqae1[i]):od:
Equation:=Array(1..node+nae):
j11 := Matrix(1 .. Nt, 1 .. Nt, storage = sparse);
J1:=Matrix(1 .. Nt, 1 .. Nt, storage = sparse);
for i from 1 to Nt do 
L:=indets(eqs[i])minus{t};
LL := [seq(ListTools:-Search(L[j], Vars),
             j = 1 .. nops(L))];
 LL2 := ListTools:-MakeUnique(LL);
if LL2[nops(LL2)] = 0 then
             LL := [seq(LL2[i1], i1 = 1 .. nops(LL2) - 1)]
         end if;
         if LL2[1]= 0 then LL:=[seq(LL2[i1],i1=2..nops(LL2))]: end:#print(LL);
varsdir:=[seq(Vars[LL[i]]=uu[LL[i]]+uu_old[LL[i]],i=1..nops(LL))];
varsold:=[seq(Vars[LL[i]]=uu_old[LL[i]],i=1..nops(LL))];
if i<=node then Equation[i]:=uu[i]-deltA*subs(varsdir,eqs[i])-deltA*0*subs(varsold,eqs[i]):
else Equation[i]:=subs(varsdir,eqs[i]):


end:
od:
Vars1 := [seq(uu[i], i = 1 .. Nt)];
    
    for i to Nt do
        L := indets(Equation[i]);
        LL := [seq(ListTools:-Search(L[j], Vars1),
            j = 1 .. nops(L))];
        LL2 := ListTools:-MakeUnique(LL);
        if LL2[nops(LL2)] = 0 then
            LL := [seq(LL2[i], i = 1 .. nops(LL2) - 1)]
        end if;
        if LL2[1]= 0 then LL:=[seq(LL2[i],i=2..nops(LL2))]: end:

        for j to nops(LL) do
            #if LL[j]>=i then j11[i, LL[j]] := diff(Eqss[i], uu[LL[j]]);end:
            #j11[i, LL[j]] := diff(Eqss[i], uu[LL[j]]);
            j11[i, LL[j]] := diff(Equation[i], uu[LL[j]]);
        end do;
    end do;

  
F := unapply(<seq(evalf(Equation[i]), i = 1 .. node + nae)>, uu,
        uu_old, deltA);#print(F);

    Jac := unapply(j11, uu, uu_old, deltA);#print(Jac);
    ic0 := Vector(node + nae,
        [seq(rhs(ICs[i]), i = 1 .. node + nae)],
        datatype = float[8]);
    y0 := Vector(Nt, datatype = float[8]);
    yold := y0;
    err := 1;
    singlestep := proc(ic0,y00, h, node, nae,iter,Anz,NumericA)
global ncons;
        local k, db, ynew, err, yold, i, dt, ic1, ydb,Nt,J1,F1;
            Nt := node + nae;
            yold:=copy(y00):
            db:=copy(y00):
            err := 1;
            
            for k to iter while 1e-14 < err do
            
            F1:=-convert(F(yold, ic0, h),Matrix,storage=sparse):
            db:=convert(LUSolve(Nt,1,Anz,NumericA,F1),Vector):
                ynew := yold + db;
                err := evalhf(eval(LinearAlgebra:-Norm(db, Nt)/LinearAlgebra:-Norm(ynew, Nt)));
                yold := ynew;
            end do;
            
            ic1 := ic0 + ynew;
            end proc;
    tol := atol;atol1:=atol:rtol1:=10*atol1:
            count:=0;
            
            
  		yold := Vector(Nt, 0, datatype = float[8]):
               
            try
            count:=count+1:
            tt[count]:=0.0:t11:=tt[count]: 
            J1:= convert(Jac(yold, ic0,0.0),Matrix,datatype=float[8]):
            (Anz,NumericA):=LUDecomp(Nt,Nt,J1):
           
            YY[count] := singlestep(ic0,yold, 0.0, node, nae,11,Anz,NumericA);
            catch:
            tcont:=0:
            end try:
		
            dt:=hinit:
            kjac:=0:
            kmod1:=kmod:
            
            while t11<tf and count<=Ntot and tcont=1 do
            try
            count:=count+1: 
            tt[count]:=tt[count-1]+dt:
            t11:=tt[count]:
            
            if kjac = 0 then 
            J1:= convert(Jac(yold, YY[count-1], dt),Matrix,datatype=float[8]):
            (Anz,NumericA):=LUDecomp(Nt,Nt,J1):end:
            #yh:=singlestep(YY[count-1],yold, dt, node, nae,iter,Anz,NumericA):#print(yh);
            if count=2 then yh:=singlestep(YY[count-1],yold, dt, node, nae,iter,Anz,NumericA):
            else
            yh:=singlestep(YY[count-1],(YY[count-1]-YY[count-2])*dt/hold, dt, node, nae,iter,Anz,NumericA):
            end:
            ymid:=singlestep(YY[count-1], (yh-YY[count-1])/2,dt/2.0, node, nae,iter,Anz,NumericA):
             yh2:=singlestep(ymid,yh-ymid, dt/2.0, node, nae,iter,Anz,NumericA):
            YY[count]:=2*yh2-yh:#print(count,YY[count]): 

            
            
            ee := Vector(node + nae, [seq(abs(yh[i] - yh2[i])
                /(atol1 + rtol1*abs(YY[count][i])),
                i = 1 .. node + nae)]);
                errL := max(atol1,
                LinearAlgebra:-Norm(ee, 2)/(node + nae));
            if errL<1.0 then 
            growth := min(3.0,max(0.3,0.9*(1/errL)^(1/(1+1))));
            hold:=dt:
            dt:=min(tf-t11,hmax,growth*dt):#print(dt):
            
            
            if errL>0.1 then kjac:=0: else kjac:=1:end:
            else dt:=dt/4: fstep:=fstep+1:#print("step failure at",t11):
            kjac:=0:
            count:=count-1:
             if count=1 then t11:=0: else t11:=tt[count-1]:end:
            if dt<hinit/1000 then tcont:=0:print("cannot integrate past",t11)end:
            end:
            catch:
            
            tcont:=0:
            end try:
            end:
            
if tt[count]=0 then print("Simulation was stopped due to initialization failure"):tcont:=1:end:
if tt[count]<tf and tcont=0 then print("Simulation was stopped before the final time due to a state becoming complex"):end:           
if tt[count]<tf and tcont=1 then print("Simulation was stopped before the final time, increase total time steps or change other solver options"):end: 
if fstep >0 and tcont=1 then print("Integration completed, but number of failed steps=",fstep):end:
            
#[[seq([tt[i1],YY[i1]],i1=1..count)],count];
solf:=Vector(count):
for i from 1 to count do solf[i]:=[t=tt[i],seq(Vars[j]=YY[i][j],j=1..Nt)]:od:
solf;

end proc;








RadIIADAE := proc(Eqode, Eqae1, ICs,solveroptions)
global ncons;
local node, nae,  j, eqs,Vars1, varsold, n1,count,
varsdir, i, ss, Equation, k, F, Vars, j11, L, LL, LL2,
Jac, ic0, y0, yold, err, dt, singlestep, tol, Yf,t11,Nt,iter1,icb,
dt1,Anz,NumericA,J1,ymid,yh,yh2,yR,tt,ee,atol1,rtol1,errL,growth,nextH,hold,kmod1,i1,tcont,YY,
tf,atol,hinit,hmax,Ntot,iter,kmod,kjac,fstep,varsdir1,solf,varsdirint;
optimize;
with(LinearAlgebra):

tf:=rhs(solveroptions[1]):atol:=rhs(solveroptions[2]):hinit:=rhs(solveroptions[3]):hmax:=rhs(solveroptions[4]):
Ntot:=rhs(solveroptions[5]):iter:=rhs(solveroptions[6]):


YY:=Array(1..Ntot+1):
    tcont:=1;
    node := nops(Eqode);
    nae := nops(Eqae1);
    Nt:=node+nae:
    Vars := [seq(lhs(ICs[i]), i = 1 .. Nt)];
   kjac:=0:
   fstep:=0: 
yh:=Vector(2*Nt,datatype=float[8]):
yh2:=Vector(2*Nt,datatype=float[8]):
ymid:=Vector(2*Nt,datatype=float[8]):
eqs:=Array(1..node+nae):
for i from 1 to node do eqs[i]:=rhs(Eqode[i]):od:
for i from 1 to nae do eqs[i+node]:=lhs(Eqae1[i]) - rhs(Eqae1[i]):od:
Equation:=Array(1..2*Nt):
j11 := Matrix(1 .. 2*Nt, 1 .. 2*Nt, storage = sparse);
J1:=Matrix(1 .. 2*Nt, 1 .. 2*Nt, storage = sparse);
for i from 1 to Nt do 
L:=indets(eqs[i])minus{t};
LL := [seq(ListTools:-Search(L[j], Vars),
             j = 1 .. nops(L))];
 LL2 := ListTools:-MakeUnique(LL);
if LL2[nops(LL2)] = 0 then
             LL := [seq(LL2[i1], i1 = 1 .. nops(LL2) - 1)]
         end if;
         if LL2[1]= 0 then LL:=[seq(LL2[i1],i1=2..nops(LL2))]: end:#print(LL);
varsdir:=[seq(Vars[LL[i]]=uu[LL[i]]+uu_old[LL[i]],i=1..nops(LL))];
varsdirint:=[seq(Vars[LL[i]]=uu[LL[i]+Nt]+uu_old[LL[i]],i=1..nops(LL))];
#varsold:=[seq(Vars[LL[i]]=uu_old[LL[i]],i=1..nops(LL))];
if i<=node then 
 Equation[i]:=5/2*uu[i]-9/2*uu[i+Nt]-deltA*subs(varsdir,eqs[i]):
 Equation[i+Nt]:=1/2*uu[i]+3/2*uu[i+Nt]-deltA*subs(varsdirint,eqs[i]):
else Equation[i]:=subs(varsdir,eqs[i]): Equation[i+Nt]:=subs(varsdirint,eqs[i]):

 #Equation[i]:=uu[i]-deltA/4*subs(varsdir,eqs[i])-deltA*3/4*subs(varsdirint,eqs[i]):
 #Equation[i+Nt]:=uu[i+Nt]+deltA/12*subs(varsdir,eqs[i])-deltA*5/12*subs(varsdirint,eqs[i]):
#else Equation[i]:=subs(varsdir,eqs[i]): Equation[i+Nt]:=subs(varsdirint,eqs[i]):


end:
od:
Vars1 := [seq(uu[i], i = 1 .. 2*Nt)];
    
    for i to 2*Nt do
        L := indets(Equation[i]);
        LL := [seq(ListTools:-Search(L[j], Vars1),
            j = 1 .. nops(L))];
        LL2 := ListTools:-MakeUnique(LL);
        if LL2[nops(LL2)] = 0 then
            LL := [seq(LL2[i], i = 1 .. nops(LL2) - 1)]
        end if;
        if LL2[1]= 0 then LL:=[seq(LL2[i],i=2..nops(LL2))]: end:

        for j to nops(LL) do
            #if LL[j]>=i then j11[i, LL[j]] := diff(Eqss[i], uu[LL[j]]);end:
            #j11[i, LL[j]] := diff(Eqss[i], uu[LL[j]]);
            j11[i, LL[j]] := diff(Equation[i], uu[LL[j]]);
        end do;
    end do;

  
F := unapply(<seq(evalf(Equation[i]), i = 1 .. 2*node + 2*nae)>, uu,
        uu_old, deltA);

    Jac := unapply(j11, uu, uu_old, deltA);
    ic0 := Vector(2*node +2* nae,
        [seq(rhs(ICs[i]), i = 1 .. node + nae),seq(rhs(ICs[i]), i = 1 .. node + nae)],
        datatype = float[8]);
    y0 := Vector(2*Nt, datatype = float[8]);
    yold := copy(y0);
    err := 1;
#print(F(yold,ic0,0.0),Jac(yold,ic0,0.0)):
    singlestep := proc(ic0,y00, h, node, nae,iter,Anz,NumericA)
global ncons;
        local k, db, ynew, err, yold, i, dt, ic1, ydb,Nt,J1,F1;
            Nt := node + nae;
            yold:=copy(y00):
            db:=copy(y00):
            err := 1;
            
            for k to iter while 1e-14 < err do
            
            F1:=-convert(F(yold, ic0, h),Matrix,storage=sparse):
            db:=convert(LUSolve(2*Nt,1,Anz,NumericA,F1),Vector):
                ynew := yold + db;
                err := evalhf(eval(LinearAlgebra:-Norm(db, 2*Nt)/LinearAlgebra:-Norm(ynew, 2*Nt)));
                yold := ynew;
            end do;
            
            ic1 := ic0 + ynew;
            end proc;
    tol := atol;atol1:=atol:rtol1:=10*atol1:
            count:=0;
            
            
  		yold := Vector(Nt*2, 0, datatype = float[8]):
               
            try
            count:=count+1:
            tt[count]:=0.0:t11:=tt[count]: 
            J1:= convert(Jac(yold, ic0,0.0),Matrix,datatype=float[8]):
            (Anz,NumericA):=LUDecomp(2*Nt,2*Nt,J1):
           
            YY[count] := singlestep(ic0,yold, 0.0, node, nae,11,Anz,NumericA);
            catch:
            tcont:=0:
            end try:
		
            dt:=hinit:
            kjac:=0:
            kmod1:=kmod:
            
            while t11<tf and count<=Ntot and tcont=1 do
            try
            count:=count+1: 
            tt[count]:=tt[count-1]+dt:
            t11:=tt[count]:
            
            if kjac = 0 then 
            J1:= convert(Jac(yold, YY[count-1], dt),Matrix,datatype=float[8]):
            (Anz,NumericA):=LUDecomp(2*Nt,2*Nt,J1):end:
            #yh:=singlestep(YY[count-1],yold, dt, node, nae,iter,Anz,NumericA):#print(yh);
            if count=2 then yh:=singlestep(YY[count-1],yold, dt, node, nae,iter,Anz,NumericA):
            else
            yh:=singlestep(YY[count-1],(YY[count-1]-YY[count-2])*dt/hold, dt, node, nae,iter,Anz,NumericA):
            end:
            ymid:=singlestep(YY[count-1], (yh-YY[count-1])/2,dt/2.0, node, nae,iter,Anz,NumericA):
             yh2:=singlestep(ymid,yh-ymid, dt/2.0, node, nae,iter,Anz,NumericA):
            YY[count]:=8/7*yh2-yh/7: 

            
            
            ee := Vector(node + nae, [seq(abs(yh[i] - yh2[i])
                /(atol1 + rtol1*abs(YY[count][i])),
                i = 1 .. node + nae)]);
                errL := max(atol1,
                1/7*LinearAlgebra:-Norm(ee, 2)/(node + nae));
            if errL<1.0 then 
            growth := min(3.0,max(0.3,0.9*(1/errL)^(1/(3+1))));
            hold:=dt:
            dt:=min(tf-t11,hmax,growth*dt):#print(dt):
            
            
            if errL>0.1 then kjac:=0: else kjac:=1:end:
            else dt:=dt/4: fstep:=fstep+1:#print("step failure at",t11):
            kjac:=0:
            count:=count-1:
             if count=1 then t11:=0: else t11:=tt[count-1]:end:
            if dt<hinit/1000 then tcont:=0:print("cannot integrate past",t11)end:
            end:
            catch:
            
            tcont:=0:
            end try:
            end:
            
if tt[count]=0 then print("Simulation was stopped due to initialization failure"):tcont:=1:end:
if tt[count]<tf and tcont=0 then print("Simulation was stopped before the final time due to a state becoming complex"):end:           
if tt[count]<tf and tcont=1 then print("Simulation was stopped before the final time, increase total time steps or change other solver options"):end: 
if fstep >0 and tcont=1 then print("Integration completed, but number of failed steps=",fstep):end:
            
#[[seq([tt[i1],YY[i1]],i1=1..count)],count];
solf:=Vector(count):
for i from 1 to count do solf[i]:=[t=tt[i],seq(Vars[j]=YY[i][j],j=1..Nt)]:od:
solf;

end proc;