(*  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;
ApproxQ     =  0;
ErrReq      =  10^(-10);  
zRintegral  = 10;

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


(************************************************************************)
(*                   Define initial guess                               *)
(************************************************************************)
gamma = 1;
temp[1] = (1-2*sigma*gamma)/gamma;
temp[2] = (1-sigma*gamma)/gamma;
u[0]  = sigma + temp[1]*Exp[-gamma*z] - temp[2]*Exp[-2*gamma*z];

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

(************************************************************************)
(*              Defines the auxiliary linear operator                   *)
(************************************************************************)
L[u_] := D[u,{z,3}] - gamma^2 * D[u,z];

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

(* Set optimal c0 and sigma  *) 
c0    =-5/4 ;
sigma = 3;
Print[" c0  =  ",c0, "    sigma  =  ",sigma];

(*  Gain up to 10th-order HAM approxiamtion *)
BVPh[1,10];