(*  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;
NgetErr     =  1;

(************************************************************************)
(*                     Define the governing equation                    *)
(************************************************************************)  
f[z_,u_,lambda_] := D[u,{z,4}] + alpha*(z*D[u,{z,3}] + 3*D[u,{z,2}])  + 
            R( u*D[u,{z,3}] - D[u,z]*D[u,{z,2}] );  
alpha  =  3/2;
R      =  -11;
 
(***********************************************************************)
(*                     Define Boundary conditions                      *)
(***********************************************************************)
zR  = 1;
OrderEQ  = 4;
BC[1,z_,u_,lambda_] := Limit[u, z->0] ; 
BC[2,z_,u_,lambda_] := Limit[D[u,{z,2}], z->0 ];
BC[3,z_,u_,lambda_] := u - 1 /.  z->zR;
BC[4,z_,u_,lambda_] := D[u,z] /.  z->zR;

(************************************************************************)
(*                   Define initial guess                               *)
(************************************************************************)
u[0]  = sigma*z +(5-4*sigma)/2*z^3-(3-2*sigma)/2*z^5;
sigma = .;

(************************************************************************)
(*                   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]];

(*  Set c0 and sigma  *)
c0    =-1/2;
sigma = 2;
Print[" c0  =  ",c0, "    sigma  =  ",sigma];

(*  Gain HAM approximations by the 3rd-order iteration approach  *)
iter[1,40,3];