[R] help about solving two equations

[Shaoqiong Zhao]

I have two matrix s1 and s2, each of them is 1000*1.
and I have two equations:
digamma(p)-digamma(p+q)=s1,
digamma(q)-digamma(p+q)=s2,
and I want to sovle these two equations to get the value of x and y, which are also two 1000*1 matrices.

I write a program like this:

f <- function(x) {
p<- x[1]; q <- x[2]; 
 ((digamma(p)-digamma(p+q)-s1[2,]) )^2 +((digamma(q)-digamma(p+q)-s2[2,]) )^2
 }
s <- 1:10/10
g <- expand.grid(p = s, q = s)
idx <- which.min(apply(g, 1, f))
idx
g[idx,]     

I am not sure if this is correct and also this can only solve one row. How to get the whole 1000 rows of p and q?

Thanks.

Annie