matlab: bvp with constant vector in function definition

Question

The problem is quite simple:

$$\frac{d^2}{dr^2}U(r)=-4\pi rn(r)$$

The boundary-conditions are $U(0)=0$ and $U(r_{max})=c$, where c is a constant. I've tried:

    f = @(r,U) [ u(2); -4*pi.*r.*n];  
    g = @(ya,yb) [ya(1); yb(1)-c];

    solinit=bvpinit(linspace(0,r_max,10),[0 0]);
    sol = bvp4c(f,g,solinit);

For $n$ being a constant it works. If $n$ is a vector I get the error:

Error using vertcat Dimensions of matrices being concatenated are not consistent.

Is there any way to use vectors in my function definition?