(*  Install BVPh 1.0 for Mathematica 5.2.                           *)
(*  For Mathematica 8.0, please replace BVPh1_0.txt by BVPh1_1.txt  *)
<<BVPh1_0.txt;  

(************************************************************************)
(*            Define the physical and control parameters                *)
(*    TypeEQ = 1    ->   BVP in a finite domain [0,a]                   *)
(*    TypeEQ = 2    ->   eigenvalue problem in a finite domain [0,a]    *)  
(*    TypeL  = 1    ->   Chebyshev polynomial as base function          *)
(*    TypeL  = 2    ->   Hybrid-base approximation                      *)
(*         TypeBase = 1 ->   sine + polynomial                          *)
(*         TypeBase = 2 ->   cosine + polynomial                        *)
(*    ApproxQ = 0   ->   do NOT approximate the function                *)
(*    Appprox = 1   ->   approximate the function                       *)
(************************************************************************)
TypeEQ      =  1;
TypeL       =  1;
TypeBase    =  2;
ApproxQ     =  0;


(************************************************************************)
(*                     Define the governing equation                    *)
(************************************************************************)  
f[z_,u_,lambda_] := D[u,{z,4}] - beta * z * (1 + u^2);
beta = 10 ;


(***********************************************************************)
(*                     Define Boundary conditions                      *)
(***********************************************************************)
zR       = 1;
OrderEQ  = 4;
alpha    = 1/5 ;
BC[1,z_,u_,lambda_] := Limit[u, z->0 ];
BC[2,z_,u_,lambda_] := D[u,z]  /.  z->zR;
BC[3,z_,u_,lambda_] := D[u,{z,2}]  /.  z->zR;
BC[4,z_,u_,lambda_] := Module[{temp},
                      temp[1] = D[u,{z,2}] /. z->0;
                      temp[2] = D[u,{z,2}] /. z-> alpha;
                      temp[1] - temp[2] //Expand
                     ];


(************************************************************************)
(*                   Define initial guess                               *)
(************************************************************************)
U[0]  = u[0];
u[0]  = sigma/(2*alpha-3)*((6*alpha-8)*z+6*(1-alpha)*z^2+2*alpha*z^3-z^4);
sigma  = .;


(************************************************************************)
(*                   Define output term                                 *)
(************************************************************************)
output[z_,u_,k_]:= Print["output           = ", N[ u[k] /. z->1, 24] ];


(************************************************************************)
(*              Defines the auxiliary linear operator                   *)
(************************************************************************)
L[f_] := Module[{temp,numA,numB,i},
If[TypeL == 1, 
   temp[1] = D[f,{z,OrderEQ}],
   numA = IntegerPart[OrderEQ/2];
   numB = OrderEQ - 2* numA//Expand; 
   temp[0] = D[f,{z,numB}];
   For[i=1, i<=numA, i++,
       temp[1] = D[temp[0],{z,2}] + (kappa*omega[i])^2*temp[0];
       temp[0] = temp[1];  
      ];
   ];
temp[1]//Expand
];


(*  Print input  and control parameters  *)
PrintInput[u[z]];

(* Gain 3rd-order HAM approximation  *)
BVPh[1,3];

(*  Gain the squared residual of governing equation at 3rd-order approx. *)
GetErr[3];

(* Gain optimal values of c0 and sigma  *)
res=Minimize[Err[3],{c0,sigma}];
Print["Minimum square residual = ", res];

(*  Set c0 and sigma according to the above result  *)
c0=-1;
sigma = 62/100;
Print["c0 = ",c0,"    sigma  =  ",sigma];

(*  Gain 10th-order HAM approximation *)
BVPh[1,10];