Question: How to rectify the error while using HPM to coupled differential equations

Dear maple user how to rectify the error  in solving the coupled differential equation using homotropy perturbation method and direct differentiations and compare the two result by plotting the graphs :

restart:
with(PDEtools):
L:=4:Nb:=1:Nt:=1:#k is some constant
f(x):=sum((p^i)*f[i](x),i=0..L):  
g(x):=sum((p^i)*g[i](x),i=0..L):
HO1:=(1-p)*(diff(f(x),x,x))+p*(((1/x)*(diff(f(x),x))+Nb*((diff(f(x),x))*diff(g(x),x))+Nt(diff(f(x),x)^2))):                                                   
expand(%):                                                                                                                                          
collect(%,p):
HO2:=(1-p)*(diff(g(x),x,x))+p*(((1/x)*(diff(g(x),x))+(Nb/Nt)*((diff(f(x),x,x))+(1/x)*diff(f(x),x)))):                                                                                        
expand(%):                                                                                                               
collect(%,p):                                                                                                              
HO2:=%:                                                                                                                 
declare(f(x),g(x),prime=x): 
for i from 0 to L+1 do equa[1][i]:=coeff(HO1,p,i)=0 end  do:                                     
for i from 0 to L+1 do equa[2][i]:=coeff(HO2,p,i)=0 end  do:           
con[1][0]:=f[0](0)=(h(x)/64),(D(f[0]))(0)=0:                                                                                                                  
con[2][0]:=g[0](0)=-(k-h(x)^2/4),(D(g[0]))(0)=0: 
for j from 1 to L do:                                                                                                    
 con[1][j]:=f[j](0)=h(x),(D(f[j]))(0)=0:                                                                           
 con[2][j]:=g[j](0)=h(x),(D(g[j]))(0)=0:                                                                         
 end do;
for i from 0 to L do;                                                                                                     
dsolve({equa[1][i],con[1][i]},f[i](x));                                                                              
f[i](x):=rhs(%);                                                                                                         
 f[i](x):=evalf(%);                                                                                                          
dsolve({equa[2][i],con[2][i]},g[i](x));                                                                            
 g[i](x):=rhs(%);                                                                                                          
g[i](x):=evalf(%);                                                                                                         
end do; 
for u from 0 to L-1 do:                                                                                                
  f[u](_z1):subs(x=_z1,f[u](x));                                                                                    
 g[u](_z1):subs(x=_z1,g[u](x));                                                                                  
 f[u+1](x):=value(simplify(f[u+1](x)));                                                                         
   f[u+1](x):=simplify(%);                                                                                               
g[u+1](x):=value(simplify(g[u+1](x)));                                                                        
g[u+1](x):=simplify(%);                                                                                                
end do:  
f(x):=evalf(simplify(sum(f[n](x),n=0..L))); 
#### direct calculations                      
#direct solve the two equations in terms of f(x) and g(x) for h(x)=e^x where Nt and Nb #are parameters and its takes some values example 1,1
restart:
with(DETools):
with(plots):
with(IntegrationTools):
Nb:=1:Nt:=1:h(x):=e^x:
Eq1 := (diff(f(x),x,x))+(((1/x)*(diff(f(x),x))+Nb*((diff(f(x),x))*diff(g(x),x))+Nt(diff(f(x),x)^2))):   
Eq2 := (diff(g(x),x,x))+(((1/x)*(diff(g(x),x))+(Nb/Nt)*((diff(f(x),x,x))+(1/x)*diff(f(x),x)))):       

Cd1 := f(0) = h(x), (D(f))(0) = 0:
dsys := {Cd1, Eq1}:
dsol := dsolve(dsys, numeric, output = operator):
#dsol(.1):
plots[odeplot](dsol, [x, diff(f(x), x$1)], 0 .. 5, color = green):
Cd2 := g(0) = h(x), (D(g))(0) = 0:
dsys := {Cd1, Cd2, Eq1, Eq2}:
dsol := dsolve(dsys, numeric, output = operator):
plots[odeplot](dsol, [x, f(x)], 0 .. 5, color = red);
plots[odeplot](dsol, [x,g(x)], 0 .. 5, color = black);