ans = 
Loop over the values of the sigma-vector - thus call "vpaintegral" for each sigma value separately.
Or alternatively make your code a function and call it from a script with a single sigma-value as input.
Further, I'd recommend first to use "integral" instead of "vpaintegral" to save computation time.
And where do you compute w1 - w6 - the variables you try to plot at the end of your code ?
Works because "vpaintegral" is called for each sigma-value separately:
syms r
f = r^2;
sigma = 0:0.1:1;
arrayfun(@(sigma)vpaintegral(f,r,0,sigma),sigma)
Does not work because "vpaintegral" is called once with the complete sigma-vector:
syms r
f = r^2;
sigma = 0:0.1:1;
vpaintegral(f,r,0,sigma)
