I am having trouble solving a nonlinear ODE with Homotopy Analysis Method using BVPh2 package on Mathematica

Question

I have the following equations after performing the similarity transformation on fluid dynamic equations.

$f'''+ff''-2f'^2-Mf'-G_rg=0\\(1/P_r)g''+g'f-Rg + \delta Exp[g/(1+\epsilon g)]=0 $

I have 5 boundary conditions and my Initial guesses are; $U[1,0]=1 - Exp[-z];$ and $p*Exp[-z];$ and these initial guesses satisfy my boundary conditions.
My auxillary linear operators are;
$L[1,u_]:= D[u,{z, 3}]-D[u,z];$ and $L[2,u_]:= D[u,{z, 2}]-u;$

After running the BVPh2 package on Mathematica, this is the error I get;

NIntegrate::precw: The precision of the argument function (-0.5 E^(E^(-1. K[2])/(1. +0.5 E^Times[<<2>>]))) is less than WorkingPrecision (25.954589770191006`). >>  

I need your help to tell me how to make this work

Note: BVPh2 is a mathematica package that is use to solve ODE with the method of Homotopy Analysis Method (HAM).

(* Modify   control parameters in BVPh if necessary *)  
ErrReq=10^-6;  
(*           Define the governing equation          *)
TypeEQ = 1;
NumEQ = 2;
f[1,z_,{f_,g_},Lambda_]:=D[f,{z,3}] + D[f,{z,2}]*f - 2*D[f,z]*D[f,z]-M*D[f,z] - Gr*g;
f[2,z_,{f_,g_},Lambda_]:=(1/Pr)*D[g,{z,2}] + D[g,z]*f - Ra*g+del*Exp[g/(1+eps*g)];

(*              Define Boundary conditions          *)
NumBC = 5;
BC[1,z_,{f_,g_}]:=f/.z -> 0;
BC[2,z_,{f_,g_}]:=(D[f, z]-1)/.z -> 0;
BC[3,z_,{f_,g_}]:=D[f, z]/.z -> infinity;
BC[4,z_,{f_,g_}]:=(g-p)/.z -> 0;
BC[5,z_,{f_,g_}]:=g/.z -> infinity;
(*                Define solution interval          *)
zL[1]=0;
zR[1]=infinity;
zL[2]=0;
zR[2]=infinity;
zRintegral[1] =10;
zRintegral[2] =10;

(*                 Define initial guess     *)
U[1,0]= 1 - Exp[-z];
U[2,0]= p*Exp[-z];

(*          Define the auxiliary linear operator    *)
L[1,u_]:= D[u,{z, 3}]-D[u,z];
L[2,u_]:= D[u,{z, 2}]-u;
(*              Define physical parameters          *)
Gr= 1.0; M=1.0; Ra= 1.0; Pr= 0.7; p=1.0; eps=0.5; del=0.5;
(*                  Print input data                *)
PrintInput[{f[z],g[z]}];

(*                   Get Optimal Co                 *)
GetOptiVar[5,{},{c0[1],c0[2]}];                     
(*    Gain 10th-order HAM approximation             *)
BVPh[1, 5];