(* Install BVPh 1.0 for Mathematica 5.2. *) (* For Mathematica 8.0, please replace BVPh1_0.txt by BVPh1_1.txt *) < 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];