(* 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 = 2; TypeL = 1; TypeBase = 2; ApproxQ = 0; ErrReq = 10^(-20); Nupdate = 10000; NtermMax = 90; ComplexQ = 1; (************************************************************************) (* Define the governing equation *) (************************************************************************) L0[u_,z_] := D[u,{z,2}]-alpha^2*u; L02[u_,z_] := L0[L0[u,z],z]; U0 = 1 - z^2; U02 = D[U0,{z,2}]; f[z_,u_,lambda_] := L02[u,z] - I*alpha*R*(U0*L0[u,z]-U02*u) + I*lambda*R*L0[u,z]; alpha = 1 ; R = 100 ; (***********************************************************************) (* Define Boundary conditions *) (***********************************************************************) zR = 1; OrderEQ = 4; BC[0,z_,u_,lambda_] := Limit[u - 1 , z->0 ]; BC[1,z_,u_,lambda_] := Limit[D[u,z], z->0 ]; BC[2,z_,u_,lambda_] := Limit[D[u,{z,3}], z->0 ]; BC[3,z_,u_,lambda_] := u /. z->zR; BC[4,z_,u_,lambda_] := D[u,z] /. z->zR; (************************************************************************) (* Define initial guess *) (************************************************************************) U[0] = u[0]; u[0] = ( 1 - z^2 )^2; (************************************************************************) (* Defines the auxiliary linear operator *) (************************************************************************) L[f_] := D[f,{z,4}]; (************************************************************************) (* Define output term *) (************************************************************************) output[z_,u_,k_]:= Print["output = ",D[u[k],{z,2}] /. z->0//N]; (* Print input and control parameters *) PrintInput[u[z]]; (* Set convergence-control parameter *) c0 = (-1+I)/2; (* Use 3rd-iteration HAM approach *) iter[1,31,3];