Hompotopy perturbation method

restart; with(LinearAlgebra);
PDEtools[declare](f(x), prime = x);
PDEtools[declare](g(x), prime = x);
N := 3;
F := sum(p^i*f[i](x), i = 0 .. N);
G := sum(p^i*g[i](x), i = 0 .. N);
FEq := (1-p)*(diff(F, x$3))+p*(diff(F, x$3)-(2*(diff(F, x$1))*(diff(F, x$1))-(5/2*F)*(diff(F, x$2))+M*(diff(F, x$1))) . A);
GEq := (1-p)*(diff(G, x$2))+p*(diff(G, x$2)+(5/2*(1/(1+R)) . Pr)*F*(diff(G, x$1)))*K(L);

M := 1;
A := (((1-W)^2.5*(1-W))*`ρf`+W*`ρt`)/`ρf`;
Pr := 6.2;
K := ((1-W)*`ρf`*Cf+W*`ρt`*Cs)/`ρf`;
L := .5;
`ρf` := 998.3;
`ρt` := 3970;
W := .2;
R := 1.0;
Cs := 765;
Cf := 4182;
coeff(FEq, p, 0);
coeff(GEq, p, 0);
for i from 0 to N do Fequ[i] := coeff(FEq, p, i) = 0 end do;
for i from 0 to N do Gequ[i] := coeff(GEq, p, i) = 0 end do;
Fcond[1][0] := f[0](0) = 0, (D(f[0]))(0) = 1, (D(f[0]))(5) = 0; for j to N do Fcond[1][j] := f[j](0) = 0, (D(f[j]))(0) = 0, (D(f[j]))(5) = 0 end do;
Gcond[0] := g[0](0) = 1, g[0](5) = 0; for j to N do Gcond[j] := g[j](0) = 0, g[j](5) = 0 end do;

for i from 0 to N do dsolve({Fequ[i], Fcond[1][i]}, f[i](x)); f[i](x) := rhs(%) end do;
for i from 0 to N do dsolve({Gcond[i], Gequ[i]}, g[i](x)); g[i](x) := rhs(%) end do;
Fa := simplify(sum(f[n](x), n = 0 .. N)); dFa := diff(Fa, x); subs(x = 2.4, dFa);
Ga := simplify(sum(g[n](x), n = 0 .. N)); dGa := diff(Ga, x); subs(x = 2.4, dGa);
plot(Ga, x = 0 .. 5);
plot(dFa, x = 0 .. 5);

ND sove solution 

Eq1 := diff(F(x), x$3)-(2*(diff(F(x), x$1))*(diff(F(x), x$1))-(5/2*F(x))*(diff(F(x), x$2))+M*(diff(F(x), x$1))) . A = 0;
Eq2 := (diff(G(x), x$2)+(5/2*(1/(1+R)))*Pr*F(x)*(diff(G(x), x$1)))*K(L) = 0;
M := 1;
Pr := 6.2;
A := (((1-W)^2.5*(1-W))*`ρf`+W*`ρt`)/`ρf`;
K := ((1-W)*`ρf`*Cf+W*`ρt`*Cs)/`ρf`;
L := 1;
`ρf` := 998.3;
`ρt` := 3970;
W := .2;
R := 2.0;
Cs := 765;
Cf := 4182;

Cd1 := F(0) = 0, (D(F))(0) = 1, (D(F))(5) = 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) = 1, G(5) = 0;
dsys := {Cd1, Cd2, Eq1, Eq2};
dsol := dsolve(dsys, numeric, output = operator);
plots[odeplot](dsol, [x, G(x)], 0 .. 5, color = green)