interface(quiet=true):  
with(ListTools):
$include "ServiceProcs.M"
$include "02FastEqns.M"

Digits:=15:
jalpha:=0.8:
outfile:="BR.VomgValp.jalp08.dat":
outfile1:="BR.VomgValp.jalp08.a.dat":


dynvars:=[E3(z),E2(z),E1(z),Eprim3(z),Eprim2(z),Eprim1(z),h3(z),h2(z),h1(z)]:
solopts:=method=bvp[middefer],'range'=0..1,'approxsoln'=Ybvp,maxmesh=8192,'abserr'=Float(1,-5),'optimize'=true,output=mesh:
outopts:=method=bvp[middefer],'range'=0..1,'approxsoln'=Ybvp,maxmesh=8192,'abserr'=Float(1,-5),'optimize'=true,output=listprocedure:

Ybvp:=ReadInitApproxFromFile("11IAc.M",[z,E3(z),E2(z),E1(z),Eprim3(z),Eprim2(z),Eprim1(z),h3(z),h2(z),h1(z),c,Eomega,p]):

l:=20.0;
#########################################
#--- Continue in Ealpha
Ealpha0:=-84.00;  dEalpha:=-2.0;
for iii from 0 to 500 do
    Ealpha:=Ealpha0+iii*dEalpha:

    Ybvp:= dsolve(BVPROBLEM,numeric,dynvars,solopts):

    AppendToCurve(outfile,Ealpha,ValueOf(Ybvp,'c'),ValueOf(Ybvp,'Eomega'),iii):
    if (irem(iii,5)=0) then 
      OutputSolutionToFile(cat("IAP.Ealp",convert(Ealpha,'string'),".dat"),Ybvp):
      Y:= dsolve(BVPROBLEM,numeric,dynvars,outopts):
      OutputReconstructedSolutionToFile(cat("sol.Ealp",convert(Ealpha,'string'),".dat"),Y,Ealpha,ValueOf(Ybvp,'Eomega'),ValueOf(Ybvp,'c')):
    end if:
end do:

l:=20:
##########################################
#--- Continue in c:
c0:=ValueOf(Ybvp,'c'); dc:=-0.001;
Ybvp:=SwapIndepParam(Ybvp,'c','Ealpha',Ealpha): Ealpha:='Ealpha':
for iii from 0 to 500 do
    c:=c0+iii*dc:

    Ybvp:= dsolve(BVPROBLEM,numeric,dynvars,solopts):

    AppendToCurve(outfile,ValueOf(Ybvp,'Ealpha'),c,ValueOf(Ybvp,'Eomega'),iii):
    if (irem(iii,5)=0) then 
      OutputSolutionToFile(cat("IAP.c",convert(c,'string'),".dat"),Ybvp):
      Y:= dsolve(BVPROBLEM,numeric,dynvars,outopts):
      OutputReconstructedSolutionToFile(cat("sol.c",convert(c,'string'),".dat"),Y,ValueOf(Ybvp,'Ealpha'),ValueOf(Ybvp,'Eomega'),c):
    end if:
end do:

###########################################
#--- Continue in Ealpha:
Ealpha0:=ValueOf(Ybvp,'Ealpha'); dEalpha:=0.5;
Ybvp:=SwapIndepParam(Ybvp,'Ealpha','c',c): c:='c':
for iii from 0 to 500 do
    Ealpha:=Ealpha0+iii*dEalpha:

    Ybvp:= dsolve(BVPROBLEM,numeric,dynvars,solopts):

    AppendToCurve(outfile,Ealpha,ValueOf(Ybvp,'c'),ValueOf(Ybvp,'Eomega'),iii):
    if (irem(iii,5)=0) then 
       OutputSolutionToFile(cat("IAP.Ealp",convert(Ealpha,'string'),".dat"),Ybvp)::
       Y:= dsolve(BVPROBLEM,numeric,dynvars,outopts):
       OutputReconstructedSolutionToFile(cat("sol.Ealp",convert(Ealpha,'string'),".dat"),Y,Ealpha,ValueOf(Ybvp,'Eomega'),ValueOf(Ybvp,'c')):
    end if:
end do:

###########################################
#--- Continue in c:
c0:=ValueOf(Ybvp,'c'); dc:=0.005;
Ybvp:=SwapIndepParam(Ybvp,'c','Ealpha',Ealpha): Ealpha:='Ealpha':
for iii from 0 to 500 do
   c:=c0+iii*dc:

    Ybvp:= dsolve(BVPROBLEM,numeric,dynvars,solopts):

    AppendToCurve(outfile,ValueOf(Ybvp,'Ealpha'),c,ValueOf(Ybvp,'Eomega'),iii):
    if (irem(iii,5)=0) then 
       OutputSolutionToFile(cat("IAP.c",convert(c,'string'),".dat"),Ybvp):
       Y:= dsolve(BVPROBLEM,numeric,dynvars,outopts):
       OutputReconstructedSolutionToFile(cat("sol.c",convert(c,'string'),".dat"),Y,ValueOf(Ybvp,'Ealpha'),ValueOf(Ybvp,'Eomega'),c):
end if:
end do:

###########################################
#--- Continue in Ealpha:
Ealpha0:=ValueOf(Ybvp,'Ealpha'); dEalpha:=-0.5;
Ybvp:=SwapIndepParam(Ybvp,'Ealpha','c',c): c:='c':
for iii from 0 to 500 do
    Ealpha:=Ealpha0+iii*dEalpha:

    Ybvp:= dsolve(BVPROBLEM,numeric,dynvars,solopts):

    AppendToCurve(outfile,Ealpha,ValueOf(Ybvp,'c'),ValueOf(Ybvp,'Eomega'),iii):
    if (irem(iii,5)=0) then 
       OutputSolutionToFile(cat("IAP.Ealp",convert(Ealpha,'string'),".dat"),Ybvp):
       Y:= dsolve(BVPROBLEM,numeric,dynvars,outopts):
       OutputReconstructedSolutionToFile(cat("sol.Ealp",convert(Ealpha,'string'),".dat"),Y,Ealpha,ValueOf(Ybvp,'Eomega'),ValueOf(Ybvp,'c')):
    end if:
end do:
