alphavals = 1:0.5:20;
for K = 1 : length(alphavals);
alpha = alphavals(K);
a = ...
b = ...
c(K) = atan2(b,a);
end
plot(alphavals, c);
Note that I adjusted to use the four-quadrant inverse tan; the formula as you wrote it would be only the two-quadrant inverse tan. You could create that if you really wanted:
c(K) = atan(a/b);
The code here shows what you asked, looping through changing the value of alpha. But it isn't nearly as efficient as it could be. Instead, you can use
alpha = 1:0.5:20;
temp = 3E-3 .* alpha.^2; %notice the .^ instead of plain ^
a = 0.1 + temp;
b = 0.2 + 0.1 .* alpha - temp;
c = atan2(b, a);
plot(alpha, c);
