To do numerical integration, you need to use the integral function. Here, it is also necessary to use arrayfun to integrate over your ‘r’ vector. I changed ‘vth’ and also the plots, so ‘p_r’ is plotted as a function of ‘r’ (and ‘z’).
This runs, without errors or warnings. I will let your determine if it is correct:
t = 0:0.1:10*pi;
r = linspace (0, 1, numel (t));
z = linspace (0, 1, numel (t));
a = 0.009;
Gamm_inf = 2000;
nu = 3;
p0 = 1;
ro = 2;
vth = @(r) (Gamm_inf)./(2*pi*nu.*r) .* (1 - exp(-(a*r/(2*nu))));
p_r = p0 + ro.*arrayfun(@(r)integral(@(r)vth(r).*vth(r)./r,1E-8,r, 'ArrayValued',1),r)-(ro*a^2)/2.*(r.^2 + 4.*z.^2);
figure
plot(r, p_r);
grid
figure
plot3(r, z ,p_r);
grid on
Experiment to get the result you want.
