@Carl Love 

thank you very much! After using this code, it no longer prompts the problem that piecewise is not defined.

However, in my code. The result I output still has NAN.
How a variable BETA composed of variables can be used correctly in maple? My code is as follows

(*---------------------------------------------------------------------
This procedure converts a piecewise *expression* into a procedure that
    1. can be differentiated or partial differentiated with D,
    2. can be converted to Fortran and its derivatives can be
    converted to Fortran.

The optional 2nd argument is a list of the names (optionally with type
specifiers) that will be used as
the parameters of the procedure, and hence the potential variables of
differentiation. Its default value is the nonconstant names that
appear in the piecewise conditions with hfloat as the type.

This procedure uses a feature new to Maple 2021.
----------------------------------------------------------------------*)
`convert/pwproc`:= proc(
    P::specfunc(piecewise),
    V::list({name, name::type}):=
        [indets(indets(P, boolean), And(name, Not(constant)))[]]::~hfloat
)
option `Author: Carl Love <carl.j.love@gmail.com> 2021-May-9`;
    subsop(
        8= hfloat, #procedure return type
        unapply(
            subsindets(
                convert(
                    if op(0,P)::symbol or nops(P)::odd then P
                    else op(0,P)(op(P), op([0,1],P))
                    fi,
                    ifelse, ':-recurse'
                ),
                specfunc(ifelse), convert, `if`
            ),
            V
        )
    )
end proc
: 

Ha:=Wa*(2*pho(x,y)-6*pho(x,y)^2+4*pho(x,y)^3)+Aa*30*pho(x,y)^2*(1-pho(x,y))^2:
Hb:=Wb*(2*pho(x,y)-6*pho(x,y)^2+4*pho(x,y)^3)+Ab*30*pho(x,y)^2*(1-pho(x,y))^2:

theta__M:=Pi/10:

P:= piecewise(
abs(beta) <= Pi/4 and theta__M <= abs(beta),
1 + eps*cos(4*beta),
abs(beta) < theta__M,
(1 + eps*cos(4*theta__M))*cos(beta)/cos(theta__M),
1 + eps*cos(4*beta)
);

PP:= convert(P, pwproc);

dPP:= D[1](PP);

Digits:= 18:
CodeGeneration:-Fortran(PP);
CodeGeneration:-Fortran(dPP);

eqA:= ( PP^2*( nab(test(x,y))[i] &t nab(ph(x,y))[i])))
              +( dPP/2 * ( diff(test(x,y),x)*diff(pho(x,y),y)
              -diff(test(x,y),y)*diff(pho(x,y),x))):
eqA:=  eval(eqA,beta=t_atan2(diff(pho(x,y),y), diff(pho(x,y),x)));

*Here, pho(x,y) has not yet been defined

Is it because of this reason that Beta cannot be judged, and thus the formula selection in piecewise cannot be determined?
I have used if...then... this problem also occurs, but when I give a beta volume, in this case like Pi/4, then it can choose from definition.

@Carl Love 

thank you, and sorry for reply delay.

i just received campus license for Maple2021, i will try it tonight and hope can send you the test results in later time or tomorrow.

@Carl Love 

Very correct. . .

@Carl Love 

I have installed the SLATEC.. but it still not works.
Now, i am trying to use the way you suggested me, but i dont know how to write the subsindets with piecewise, Could you plesae send me a Example of reference? thank you

@Carl Love 

Thank you!  i will try it tomorrow. Tonight i just find a package named PPPACK (Piecewise Polynomial Package) . I dont know if it useful and i am trying to find a way to insatll it.

https://people.math.sc.edu/Burkardt/f_src/pppack/pppack.html

Anyway, i will report here when i get new result.

@Carl Love 

Many thanks for the answer !!!!    and sorry for the misspell....

I have tried above code and it is effectively helped me exclued undefeined in piecewise!!!!!!!!!!

However, I still have not been able to confirm the final result yet.

The compilation on maple is passed, but when i do the Make on g77, it prompts undefined reference to 'piecewise_'. It may be that my SLATEC does not include some required packages. I am now trying to install the environment....to determine the cause of the problem.

@Carl Love 

Many thanks for the answer and recommendation.

diffW2:= combine(convert(diff(W^2, beta), piecewise, beta));

seems to be helpful in solving this problem, however,  theta__M is a constant like(near to) Pi/6

so it cannot be assumed.

Now there are the following test results. When the code is as follows.
 

Ha:=Wa*(2*pho(x,y)-6*pho(x,y)^2+4*pho(x,y)^3)+Aa*30*pho(x,y)^2*(1-pho(x,y))^2:
Hb:=Wb*(2*pho(x,y)-6*pho(x,y)^2+4*pho(x,y)^3)+Ab*30*pho(x,y)^2*(1-pho(x,y))^2:

theta__M :=Pi/6:
W:= piecewise(
       beta < Pi/2-theta__M,
            0,
       beta <= Pi/2+theta__M, 
            (1+eps*cos(4*(theta__M)))*cos(beta)/cos(theta__M),
       beta <= Pi-theta__M, 
            1+eps*cos(4*beta));
#W(arctan(diff(pho(x,y),y), diff(pho(x,y),x)));

noise :=alpha*tr(el)*16*pho(x,y)^2*(1-pho(x,y))^2*((1-co(x,y))*Ha+co(x,y)*Hb):

diffW2 :=Combin(convert(diff(W^2,beta),piecewise,beta));

feqB:=  1/dt*ElementInt(test(x,y)*(ph(x,y)-pho(x,y)))
     = 
-M*ElementInt( test(x,y)*((1-co(x,y))*Ha+co(x,y)*Hb))
     +M*ElementInt( test(x,y)*noise)
+M*gamma*(-ElementInt( W^2*( nab(test(x,y))[i] &t nab(ph(x,y))[i])))
              +ElementInt( diffW2/2 * ( diff(test(x,y),x)*diff(pho(x,y),y)
              -diff(test(x,y),y)*diff(pho(x,y),x))):

lfeq:=  subs(beta=(t_atan2(diff(pho(x,y),y), diff(pho(x,y),x))),feqB);

The error prompted is as follows

Error, (in PiecewiseTools:-Is) Wrong kind of parameters in piecewise

After making the following changes to W
 

Ha:=Wa*(2*pho(x,y)-6*pho(x,y)^2+4*pho(x,y)^3)+Aa*30*pho(x,y)^2*(1-pho(x,y))^2:
Hb:=Wb*(2*pho(x,y)-6*pho(x,y)^2+4*pho(x,y)^3)+Ab*30*pho(x,y)^2*(1-pho(x,y))^2:

theta__M :=Pi/6:
W:= piecewise(
       
            
       beta <= Pi/2+theta__M, 
            (1+eps*cos(4*(theta__M)))*cos(beta)/cos(theta__M),
       beta <= Pi-theta__M, 
            1+eps*cos(4*beta));
#W(arctan(diff(pho(x,y),y), diff(pho(x,y),x)));

noise :=alpha*tr(el)*16*pho(x,y)^2*(1-pho(x,y))^2*((1-co(x,y))*Ha+co(x,y)*Hb):

diffW2 :=Combin(convert(diff(W^2,beta),piecewise,beta));

feqB:=  1/dt*ElementInt(test(x,y)*(ph(x,y)-pho(x,y)))
     = 
-M*ElementInt( test(x,y)*((1-co(x,y))*Ha+co(x,y)*Hb))
     +M*ElementInt( test(x,y)*noise)
+M*gamma*(-ElementInt( W^2*( nab(test(x,y))[i] &t nab(ph(x,y))[i])))
              +ElementInt( diffW2/2 * ( diff(test(x,y),x)*diff(pho(x,y),y)
              -diff(test(x,y),y)*diff(pho(x,y),x))):

lfeq:=  subs(beta=(t_atan2(diff(pho(x,y),y), diff(pho(x,y),x))),feqB);

There is no more error message, but undefined and 0 appear in diffw2

diffW2 := Combin(piecewise(beta < (2*Pi)/3, -(8*(1 - eps/2)^2*cos(beta)*sin(beta))/3, beta = (2*Pi)/3, undefined, beta < (5*Pi)/6, -8*(1 + eps*cos(4*beta))*eps*sin(4*beta), beta = (5*Pi)/6, undefined, (5*Pi)/6 < beta, 0))

Remarks, ElementInt is finite element FEM formula.

reference and citation：https://www.mech.kth.se/~gustava/femLego/test.html

@Carl Love 

Thank you for the answer, The description I provided earlier is not complete, W will also be taken into other calculations. code show as below

feqB:=  1/dt*ElementInt(test(x,y)*(ph(x,y)-pho(x,y)))
     = 
-M*ElementInt( test(x,y)*((1-co(x,y))*Ha+co(x,y)*Hb))
     +M*ElementInt( test(x,y)*noise)
+M*gamma*(-ElementInt( W^2*( nab(test(x,y))[i] &t nab(ph(x,y))[i])))
              +ElementInt( diff(W^2,beta)/2 * ( diff(test(x,y),x)*diff(pho(x,y),y)
              -diff(test(x,y),y)*diff(pho(x,y),x))):
lfeq:=  subs(beta=t_atan2(diff(pho(x,y),y), diff(pho(x,y),x)),feqB);

When W is improved according to your writing, the maple compilation problem is gone.

but the following problem occurred.

+ElementInt( diff(W^2,beta)/2 * ( diff(test(x,y),x)*diff(pho(x,y),y)
              -diff(test(x,y),y)*diff(pho(x,y),x))):
 

become to  ElementInt(0)

In the initial W=1+eps*cos(4*beta), the result of this part is,

ElementInt(-4*(1 + eps*cos(4*t_atan2(diff(pho(x, y), y), diff(pho(x, y), x))))*eps*sin(4*t_atan2(diff(pho(x, y), y), diff(pho(x, y), x)))*(diff(test(x, y), x)*diff(pho(x, y), y) - diff(test(x, y), y)*diff(pho(x, y), x)))