Is it possible to differentiate an equation with respect to a variable in VUHARD subroutine?

I am using fortran to write a VUHARD subroutine to implement an constitutive equation (Johnson-Cook model) given by (A + B*(strain**xn))*(one - T_star**xm)* (one + C*edot_log) in ABAQUS. I am getting the value of yield for each iteration by putting in the values of strain, T_star and edot_log from previous iteration. However, my code also involves three differentiation
 1). yield(km) with respect to strain
 2). yield(km) with respect to strainrate
 3). yield(km) with respect to T
 Currently, I am finding out the derivative by hand and then writing down the equations:
         dyieldDeqps(km,1) = (B*xn*(strain**(xn-one)))*

1               (one - T_star**xm) * (one + C*edot_log)

        dyieldDeqps(km,2) = (A + B*(strain**xn))*

1               (one - T_star**xm) * (C/strainrate)

        dyieldDtemp(km) = (-xm)*(one/(T - T_trans))*

1               (T_star**(xm))*(A + B*(strain**xn))

2               *(one + C*edot_log)
 This works fine, but, is there a way where I can find the derivatives directly from the yield equation?(A + B*(strain**xn))*(one - T_star**xm)* (one + C*edot_log)

 Any help or suggestion would be most helpful. I don't know how to do this?
 The reason I would like to do this is because I have another long constitutive equation that I want implement as a VUHARD. However, the derivatives of the yield equation becomes very long to write in the form given in the example. Is there a way to automatically do the procedure?
 The relevant parts of the FORTRAN subroutine:
 
     T_star = (T - T_trans)/(T_melt - T_trans)
    strainrate = eqpsRate(km)/epsdot0
    edot_log = LOG( strainrate )


    yield(km) = (A + B*(strain**xn))*(one - T_star**xm)

1               * (one + C*edot_log)



    IF( strain .le. zero) THEN
        dyieldDeqps(km,1) = zero

    ELSE

        dyieldDeqps(km,1) = (B*xn*(strain**(xn-one)))*

1               (one - T_star**xm) * (one + C*edot_log)

    ENDIF

    IF(strainrate .ge. epsdot0) then

        dyieldDeqps(km,2) = (A + B*(strain**xn))*

1               (one - T_star**xm) * (C/strainrate)

    ELSE

        dyieldDeqps(km,2) = one

    ENDIF

    IF (T .gt. T_trans) THEN

        dyieldDtemp(km) = (-xm)*(one/(T - T_trans))*