You can easily use the optim() function in R to find out the ML estimates of the parameters, indirectly. This means that you have to define the likelihood function properly and introduce some initial values, then use the optim() function to extract MLE numerically. You also define, directly, the derivatives of the likelihood function with respect to parameters and solve the system. In this way, you should compute digamma function numerically by some available methods.
Briefly, I recommend to use optim() function in R to catch your estimates.