(*  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);  
NgetErr     =  100;
zRintegral  =  10;


(************************************************************************)
(*                     Define the governing equation                    *)
(************************************************************************) 
f[z_,u_,lambda_] := Module[{temp},
temp[1] = D[u,{z,3}] + (1-t)*z/2*D[u,{z,2}] - t*(1-t)*D[u,z,t];
temp[2] = u*D[u,{z,2}] - D[u,z]^2;
temp[1] + t* temp[2] // Expand
];
 
(***********************************************************************)
(*                     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                               *)
(************************************************************************)
u[0]  =  1 - Exp[-z];

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

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


(* Define Getdelta[k]  *)
Getdelta[k_]:=Module[{temp,i},
uz[k]   = D[u[k],z]//Expand;
uzz[k]  = D[uz[k],z]//Expand;
uzzz[k] = D[uzz[k],z]//Expand;
uzuz[k] = Sum[uz[i]*uz[k-i],{i,0,k}]//Expand;
uzzu[k] = Sum[uzz[i]*u[k-i],{i,0,k}]//Expand;
uzt[k]  = D[uz[k],t]//Expand;
temp[1] = uzzz[k] + (1-t)*z/2*uzz[k]-t*(1-t)*uzt[k]//Expand;
temp[2] = t*(uzzu[k] - uzuz[k]);
delta[k] = temp[1] + temp[2]//Expand;
];


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

(*  Set convergence-control parameter c0 *)
c0  = -1/4;
Print["c0  =  ",c0];

(*  Gain 6th-order HAM approximation *)
Print[" c0  =  ",c0];
BVPh[1,10];


(* Calculate the squared residual  *)
For[k=2, k<=10, k=k+2,
GetErr[k];
err[k] = Integrate[Err[k],{t,0,1}];
Print[" k  =  ", k, "   Squared residual =  ", err[k]//N];
];