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


(************************************************************************)
(*                     Define the governing equation                    *)
(************************************************************************) 
L0[z_, u_] := Sqrt[1+z^2]*D[u,{z,2}] + Cos[Pi*z]*D[u,z]/z;  
F[z_,u_]  := Exp[u]/(1+z^2) + (1+z)*Sin[u];
g[z_]     := Sin[z^2 + Exp[-z]]; 
f[z_,u_,lambda_]:= L0[z,u] + lambda*F[z,u] - g[z]; 

 
(***********************************************************************)
(*                     Define Boundary conditions                      *)
(***********************************************************************)
zR       = Pi;
OrderEQ  = 2;
BC[0,z_,u_,lambda_] := Limit[u - A, z->0] ; 
BC[1,z_,u_,lambda_] := Limit[D[u,z], z->0 ];
BC[2,z_,u_,lambda_] := u - D[u,z] - 3/5 /.  z-> zR;
A   =  1/2;

(************************************************************************)
(*                   Define initial guess                               *)
(************************************************************************)
u[0]  = (5*A+3)/10 + (5*A-3)/10*Cos[kappa*z]; 
kappa = 1;


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

(************************************************************************)
(*              Defines the auxiliary linear operator                   *)
(************************************************************************)
omega[1] = 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 approach  *)
iter[1,6,3];