I'm not sure if I've interpreted your requirements correctly, but what about the following:
%% Inputs
IR = 1e-3; %initial radius [m]
Rho = 1000; %density [kg/m^3]
T = 23.5; %temperature [celsius]
% Loop 1
aInitial = 1;
aStep = 1;
aMax = 25;
VMin = 1e-12; %minimum volume [m^3]
VMax = 2e-9; %maximum volume [m^3]
% Loop 2
bInitial = 1;
bStep = 1;
bMax = 25;
ReHuMin = 0.0; %minimum relative humidity [m]
ReHuMax = 0.9; %maximum relative humidity [m]
% Loop 3
c = 1;
tInitial = 0; %initial time [s]
tStep = 1; %final time [s]
D_T = 2.5e-4*exp(-684.15/(T+273.15)); %coefficient [m^2/s]
c_sat = (9.99e-7)*T^3 - (6.94e-5)*T^2 + (3.2e-3)*T - 2.87e-2; %concentration [kg/m^3]
%% Computing
for a = aInitial:aStep:aMax
Va = VMin + (VMax - VMin)*(a-aInitial)/(aMax);
for b = bInitial:bStep:bMax
V(a,b) = Va; %Initial volume
Vol(a,b) = Va;
h(a,b) = ((sqrt(pi^2*IR^6 + 9*(Vol(a,b))^2) + 3*(Vol(a,b)))^(2/3) - pi^(2/3)*IR^2)/(pi^(1/3)*(sqrt(pi^2*IR^6 + 9*(Vol(a,b))^2) + 3*(Vol(a,b)))^(1/3));
CAR(a,b) = 2*atan(h(a,b)/IR);
CAD(a,b) = CAR(a,b)*180/pi;
RH(a,b) = ReHuMin + (ReHuMax - ReHuMin)*(b - bInitial)/(bMax);
timeSec(a,b) = tInitial;
Vold = 0; c = 0;
while Vol(a,b)>VMin && c<10000
Vold = Vol(a,b);
M_dot(a,b) = -pi*IR*D_T*(1 - RH(a,b))*c_sat*(0.27*CAR(a,b)^2+1.30); %mass flow
Mkg(a,b) = M_dot(a,b)*tStep; %mass loss at each time step [kg]
Vm3(a,b) = Mkg(a,b)/Rho; %volume loss [m^3]
Vol(a,b) = Vol(a,b) + Vm3(a,b); %new volume [m^3]
Vol(a,b) = max(Vol(a,b), VMin);
h(a,b) = ((sqrt(pi^2*IR^6 + 9*(Vol(a,b))^2) + 3*(Vol(a,b)))^(2/3) - pi^(2/3)*IR^2)/(pi^(1/3)*(sqrt(pi^2*IR^6 + 9*(Vol(a,b))^2) + 3*(Vol(a,b)))^(1/3)); %new height [m]
CAR(a,b) = 2*atan(h(a,b)/IR); %new angle [radians]
CAD(a,b) = CAR(a,b)*180/pi; %new angle [degrees]
timeSec(a,b) = tStep*c;
c = c + 1;
end
EvaporationTime(a,b) = timeSec(a,b);
end
end
%% Plotting
contourf(RH*100,V*1e9,EvaporationTime) % plot against Initial volume not final volume
k = colorbar;
k.Label.String = 'Evaporation Time (s)';
xlabel('Relative Humidity (%)')
ylabel('Initial Volume (mm^3)')
figure
surf(RH*100,V*1e9,EvaporationTime)
k = colorbar;
k.Label.String = 'Evaporation Time (s)';
xlabel('Relative Humidity (%)')
ylabel('Initial Volume (mm^3)')
zlabel('Evaporation Time (s)')
