Using glmfit for logistic regression, how to do it?

I have a set of data points (r, x1, x2, x3, D) where r, x1, x2, x3 describe the location of the point in 3D space and D is the difficulty associated with each point. I am trying to create a series of 3D curved surface to represent each difficulty D by doing a logistic regression to yield a unique model of the form r = beta(1)x1 + beta(2)x2 + beta(3)x3 for each D.

I am trying to use a binomial distribution, using 0's to represent points that do not contain the D of interest, and a 1 for those that do. So for each D I am computing:

beta = glmfit([x1 x2 x3], [isD], 'binomial', 'constant', 'off');

But this does not account for the fact that I need an equation to calculate r, not a probability density function which is what I think I'm getting.

I know that a linear regression won't work because I've tried it and my data is not linear in nature. How do I go about doing a logistic regression properly?