I'm trying to solve nonlinear system of equations via homotopy perturbation method. Please help with ideas to develop the code.

restart;

de1 := (1 - p)*(diff(x(t), t $ 1) - r*x) + p*(diff(x(t), t $ 1) - r*(1 - (x + y)/K) + al*x*v);

de2 := (1 - p)*(diff(y(t), t $ 1) + m*y+u*y) + p*(diff(y(t), t $ 1) - al*x*v + m*y+u*y);

de3 := (1 - p)*(diff(v(t), t $ 1) + eta*v) + p*(diff(v(t), t $ 1) - B*y + eta*v);

ibvc := x(0) = 0.005, y(0) = 0.0007;

v(0) = 2;

r:=0.005; al:=0.002; m:=0.0001; K:=0.0138; B:=0.2; eta:=0.006; u:=0.0001;

n := 4;

X := unapply(add(g[k](t)*p^k, k = 0 .. n), x);

Y := unapply(add(h[k](t)*p^k, k = 0 .. n), x);

V := unapply(add(i[k](t)*p^k, k = 0 .. n), x);

DE1 := series(eval(de1, x = X), p = 0, n + 1);

DE2 := series(eval(de2, y = Y), p = 0, n + 1);

DE3 := series(eval(de3, v = V), p = 0, n + 1);

C1 := map(coeffs, eval([ibvc], x = X), p);

C2 := map(coeffs, eval([ibvc], y = Y), p);

C3 := map(coeffs, eval([ibvc], v = V), p);

for k from 0 to n do

    IBVC1 := select(has, C1, g[k]);

    IBVC2 := select(has, C2, h[k]);

    IBVC3 := select(has, C3, i[k]);

    s1 := dsolve({coeff(DE1, p, k), op(IBVC1)});

    s2 := dsolve({coeff(DE2, p, k), op(IBVC2)});

    s3 := dsolve({coeff(DE3, p, k), op(IBVC3)});

    h[k] := unapply(rhs(s1), t);

h[k] := unapply(rhs(s2), t);

i[k] := unapply(rhs(s3), t);

end do;

'X(t)' = X(t) + O(p^(n + 1));

'Y(t)' = Y(t) + O(p^(n + 1));

'V(t)' = V(t) + O(p^(n + 1));

How to make changes in the code?