# creates an array of "numpoints" abscisas 
# starting at "x_start" and ending at "x_end" 
#
xValues:=proc(x_start,x_end,numpoints)
local dx, i:
    dx:=(x_end-x_start)/numpoints:
    [seq(x_start+i*dx, i=0..numpoints)]:
end:


#----------------
xyValues:=proc(f::procedure,x_start,x_end,numpoints)
local dx, i:
    dx:=(x_end-x_start)/numpoints:
    Array([seq([x_start+i*dx,f(x_start+i*dx)], i=0..numpoints)]):
end:

#----------------
xyyValues:=proc(f::procedure,g::procedure,x_start,x_end,numpoints)
local dx, i:
    dx:=(x_end-x_start)/numpoints:
    Array([seq([x_start+i*dx,f(x_start+i*dx),g(x_start+i*dx)], i=0..numpoints)]):
end:

#----------------
xyyyValues:=proc(f::procedure,g::procedure,w::procedure,x_start,x_end,numpoints)
local dx, i:
    dx:=(x_end-x_start)/numpoints:
    Array([seq([x_start+i*dx,f(x_start+i*dx),g(x_start+i*dx),w(x_start+i*dx)], i=0..numpoints)]):
end:


#----------------
y02Values:=proc(f::procedure,g::procedure,w::procedure,x_start,x_end,numpoints)
local dx, i:
    dx:=(x_end-x_start)/numpoints:
    Array([seq([f(x_start+i*dx),g(x_start+i*dx),w(x_start+i*dx)], i=0..numpoints)]):
end:

#----------------
xy03Values := proc(f::procedure,g::procedure,E::procedure,n::procedure, x_start, x_end, numpoints)
local dx, i:
    dx := (x_end - x_start)/numpoints: 
    Array([seq([x_start + i*dx, f(E(x_start + i*dx),n(x_start + i*dx)),g(x_start + i*dx)], i = 0 .. numpoints)]):
end proc:


#===========================================================================


ReadInitApproxFromFile:=proc(file::string,vars::list)
# Reads data from file and prepares an initial approximation for "dsolve,bvp"
  local B,len,InitApprox:
  len:=nops(vars);
  B:=readdata(file,len):
  B:=Transpose(B):
  InitApprox:=[ seq (FlattenOnce([vars[i],B[i]]), i=1..len)]:
InitApprox;
end proc:

#-----------------
IndexOf:=proc(SOL::Matrix,Var)
# Returns the index of the component "Var" in the solution matrix "SOL" of the BVP.
  local B,wid,len,jjj,Indx:
 B:=convert(SOL[2,1],'listlist'):   wid:=nops(B[1]);
 B:=Transpose(B):                   len:=nops(B[1]);
       for jjj from 1 to wid do
         if SOL[1,1][jjj] = Var then
            Indx:=jjj;
         end if;
       end do:
Indx;
end proc:

#------------
ValueOf:=proc(SOL::Matrix,Var::symbol)
# Returns the value of a _parameter_ "Var" from the solution matix "SOL" of the BVP.
  local B,wid,len,Indx:
 B:=convert(SOL[2,1],'listlist'):   wid:=nops(B[1]);
 B:=Transpose(B):                   len:=nops(B[1]);
 Indx:=IndexOf(SOL,Var):
B[Indx,len];
end proc:


#--------------
SwapIndepParam:=proc(SOL::Matrix,oldpar::symbol,newpar::symbol,newparValue)
# Prepares an initial approximation with respect to a new continuation parameter
  local B,wid,len,Indx,oldparValue,vars,NEWSOL:
  B:=convert(SOL[2,1],'listlist'):   wid:=nops(B[1]);
  B:=Transpose(B):                   len:=nops(B[1]);
  Indx:=IndexOf(SOL,oldpar):
  vars:=SOL[1,1]:
  vars:=[seq(vars[i], i=1..Indx-1), newpar, seq (vars[i], i=Indx+1..wid)]:
  B[Indx]:=[seq(newparValue,i=1..len)]:
  NEWSOL:=[ seq (FlattenOnce([vars[i],B[i]]), i=1..wid)]:
NEWSOL;
end proc:

#-----------------
AppendToCurve:=proc(file::string,par1,par2,par3,par4)
# Does what it does. "All I wanna do is do it", Kevin the Teenager.
  local fd;
   printf("%2.5E %2.20E %2.20E %a \n",par1,par2,par3,par4);
   fd := fopen(file, APPEND):
   fprintf(fd, "%2.20E %2.20E %2.20E\n",par1,par2,par3):
   fclose(fd):
end proc:

#-----------------
OutputSolutionToFile:=proc(file::string,SOL::Matrix)
# Ouput all components of the solution matrix to file so that they it can be 
#   used as an initial approxinamation for onother "dsolve,numeric,bvp".
  local B,wid,C,len,Ba,fd: 
  B:=convert(SOL[2,1],'listlist'):   wid:=nops(B[1]);
  C:=Transpose(B):                   len:=nops(C[1]);
  Ba:=Array([seq( FlattenOnce([B[i]]), i=1..len) ]):
   fd := fopen(file, WRITE):
   fprintf(fd, "%+2.15E \n", Ba):
   fclose(fd):
end proc:

#-----------------
OutputReconstructedSolutionToFile:=proc(file::string,SOL::list,par1,par2,par3)
# Ouput all components of the solution matrix to file so that they it can be 
#   used as an initial approxinamation for onother "dsolve,numeric,bvp".
  local EA,EB,EC,HA,HB,HC,PA,PAA,XX,XB,X,EE,HH,Ealp,Eomg,cval,fd:

EA := eval(E3(z),SOL): EB := eval(E2(z),SOL): EC := eval(E1(z),SOL):
HA := eval(h3(z),SOL): HB := eval(h2(z),SOL): HC := eval(h1(z),SOL):
PA:= eval(p,Y): PAA:=PA(1):
XX:=xValues(0,1,200): XB:=xValues(1,0,200):
X:=Flatten([map(x->-(l*x)-l/PAA,XB),map(x->-l*x/PAA,XB),map(x->x*l,XX)]):
EE:=Flatten([map(x->EC(x),XB),map(x->EB(x),XB), map(x->EA(x),XX)]):
HH:=Flatten([map(x->HC(x),XB),map(x->HB(x),XB), map(x->HA(x),XX)]):
Ealp:=Flatten([map(x->par1,XB),map(x->par1,XB), map(x->par1,XX)]):
Eomg:=Flatten([map(x->par2,XB),map(x->par2,XB), map(x->par2,XX)]):
cval:=Flatten([map(x->par3,XB),map(x->par3,XB), map(x->par3,XX)]):
fd := fopen(file, WRITE):
fprintf(fd,"%+2.15E \n", Array(Transpose([X,EE,HH,Ealp,Eomg,cval]))): 
fclose(fd):
end proc:

