It looks like your error message comes from an ODE solver. It's the ODE solver that calls your function, so it's unclear what you mean by "I tried simply calling..."
Anyway, if you get this error, the problem is not with the way the function is called, it's with the way the function constructs its outputs. It's currently constructing the outputs as row vector, whereas the calling code expects column vector.
The simplest way is to preallocate the outputs as column vectors. Outputs should be preallocated anyway and with preallocation you can define the shape. So, before the loop:
steps = PI:PI+deltaq; %steps of the for loop
nsteps = numel(steps); %number of steps
V = zeros(nsteps, 1); %preallocate V rather than growing it in the loop
Dvdq = zeros(nsteps, 1); %same
F = zeros(nsteps, 1);
dFdq = zeros(nsteps, 1);
dPdq = zeros(nsteps, 1); %and preallocate output as a COLUMN vector.
for q = steps %your loop
%... rest of the code unchanged
end
Also note that using global variables is a bad idea. See Parameterizing Functions for better methods.
