#Code written by Venkat R. Subramanian and MAPLE group members, UT Austin, last updated on 11/30/2022. Released under MIT license #To run this code users need Intel's MKL libraries (Pardiso) pFactor := define_external("hw_PardisoFactor", ':-MAPLE', ':-LIB' = "linalg"): pSolve := define_external("hw_PardisoSolve", ':-MAPLE', ':-LIB' = "linalg"): pCSF := define_external("hw_SpCompressedSparseForm", ':-MAPLE', ':-LIB' = "linalg"): IMPDAEP := 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, CB,R,X,handle,db; 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(, 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]); db := Vector(Nt, datatype = float[8]); yold := y0; err := 1; singlestep := proc(ic0,y00, db,h, node, nae,iter,handle) local k, 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:=-F(yold, ic0, h): pSolve(F1,db,handle): #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]): J1:=convert(Jac(yold, ic0, 0.0),Matrix,storage=sparse,datatype=float[8]): CB,R,X:=pCSF(J1,"r",1,1,false,true,false): handle := pFactor(CB, R, X, 11): #(Anz,NumericA):=LUDecomp(Nt,Nt,J1): YY[count] := singlestep(ic0,yold, db,0.0, node, nae,11,handle);#print(YY[count],t11): catch: tcont:=0: end try: dt:=hinit: kjac:=0: kmod1:=kmod: while t110.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 dt0 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; CNDAEP := 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, CB,R,X,handle,db; 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))]; 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/2*subs(varsdir,eqs[i])-deltA/2*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(, 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]); db := Vector(Nt, datatype = float[8]); yold := y0; err := 1; singlestep := proc(ic0,y00, db,h, node, nae,iter,handle) local k, 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:=-F(yold, ic0, h): pSolve(F1,db,handle): #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]): J1:=convert(Jac(yold, ic0, 0.0),Matrix,storage=sparse,datatype=float[8]): CB,R,X:=pCSF(J1,"r",1,1,false,true,false): handle := pFactor(CB, R, X, 11): #(Anz,NumericA):=LUDecomp(Nt,Nt,J1): YY[count] := singlestep(ic0,yold, db,0.0, node, nae,11,handle);#print(YY[count],t11): catch: tcont:=0: end try: dt:=hinit: kjac:=0: kmod1:=kmod: while t110.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 dt0 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; EBDAEP := 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, CB,R,X,handle,db; 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))]; 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(, 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]); db := Vector(Nt, datatype = float[8]); yold := y0; err := 1; singlestep := proc(ic0,y00, db,h, node, nae,iter,handle) local k, 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:=-F(yold, ic0, h): pSolve(F1,db,handle): #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]): J1:=convert(Jac(yold, ic0, 0.0),Matrix,storage=sparse,datatype=float[8]): CB,R,X:=pCSF(J1,"r",1,1,false,true,false): handle := pFactor(CB, R, X, 11): #(Anz,NumericA):=LUDecomp(Nt,Nt,J1): YY[count] := singlestep(ic0,yold, db,0.0, node, nae,11,handle);#print(YY[count],t11): catch: tcont:=0: end try: dt:=hinit: kjac:=0: kmod1:=kmod: while t110.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: #t11:=tt[count-1]: if count=1 then t11:=0: else t11:=tt[count-1]:end: if dt0 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; RadIIADAEP := 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,db,CB,X,handle,R; 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(, 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]); db := 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, db,h, node, nae,iter,handle) local k, ynew, err, yold, i, dt, ic1, ydb,Nt,J1,F1; Nt := node + nae; yold:=copy(y00): err := 1; for k to iter while 1e-14 < err do F1:=-F(yold, ic0, h): pSolve(F1,db,handle): yold := yold + db; err := evalhf(eval(LinearAlgebra:-Norm(db, Nt)/LinearAlgebra:-Norm(yold, Nt))); end do; ic1 := ic0 + yold; 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,storage=sparse,datatype=float[8]): CB,R,X:=pCSF(J1,"r",1,1,false,true,false): handle := pFactor(CB, R, X, 11): #J1:= convert(Jac(yold, ic0,0.0),Matrix,datatype=float[8]): #(Anz,NumericA):=LUDecomp(2*Nt,2*Nt,J1): YY[count] := singlestep(ic0,yold, db,0.0, node, nae,11,handle); catch: tcont:=0: end try: dt:=hinit: kjac:=0: kmod1:=kmod: while t110.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 dt0 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;