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

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

Maple

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);

```
﻿