I see a couple of problems in the code:
First, ampu(n+1) means a scaler, but sqrt(au.^2+bu.^2) is a vector. So, you can form ampu as a matrix whose rows corresponds to each harmonics, like: ampu(n,:)=sqrt(au.^2+bu.^2);
Second, in plot(x,ampu), x and ampu are not in same lengths. Look, you defined x=0:5:200; which has 41 elements, but ampu has 201 elements. Make them of same lenghts. For example define x=linspace(0,200,201); or as per your requirement.
So, run the following modified code, which produces animation-like plots of ampu for 40 times. You can increase the pause to make it slower. However, if you want to see all 40 harmonics at a time then plot them outside the for loop.
t=0:0.005:1; x=linspace(0,200,201);
cu=0.5; au=zeros(1,201); bu=au; ampu=zeros(1,201); ampu(1)=0.5;
for n=1:40
au=au+1/(n*pi)*(cos(n*pi/4)-cos(n*2*pi/4)+cos(n*3*pi/4)-cos(n*4*pi/4)+cos(5*n*pi/4)-
cos(n*6*pi/4)+cos(n*7*pi/4)-cos(n*8*pi/4)) *sin(n*2*pi*t);
bu=bu+1/(n*pi)*(sin(n*2*pi/4)-sin(n*pi/4)+sin(n*4*pi/4)-sin(n*3*pi/4)+sin(n*6*pi/4)-
sin(5*n*pi/4)+sin(n*8*pi/4)-sin(n*7*pi/4)) *cos(n*2*pi*t);
ampu(n,:)=sqrt(au.^2+bu.^2);
subplot(1,2,1),plot(t,cu+au+bu),grid
subplot(1,2,2),plot(x,ampu(n,:)),grid
pause(0.01);
end
Hope this helps.
