(*  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      =  2;
TypeL       =  1;
TypeBase    =  2;
ApproxQ     =  0;
Ntruncated  =  20;


(************************************************************************)
(*                     Define the governing equation                    *)
(************************************************************************) 
f[z_,u_,lambda_] := lambda*D[u,{z,4}]-beta*z-beta*z*u^2*lambda^2;
beta = 10 ;
 
(***********************************************************************)
(*                     Define Boundary conditions                      *)
(***********************************************************************)
zR       = 1;
OrderEQ  = 4;
BC[0,z_,u_,lambda_] := u - 1  /. z -> 1  ;
BC[1,z_,u_,lambda_] := Limit[u,  z -> 0 ];
BC[2,z_,u_,lambda_] := D[u,z] /. z -> 1 ;
BC[3,z_,u_,lambda_] := D[u,{z,2}] /.  z -> 1;
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
                              ];
alpha = 1/5;

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


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


(************************************************************************)
(*              Defines the auxiliary linear operator                   *)
(************************************************************************)
omega[1] = Pi/zR;
omega[2] = Pi/zR;
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]];

(*  Use 3rd-order HAM iteration *)
iter[1,6,3];