a=1/16;c=pi^2;
t=0:.05:.5
z=zeros(size(t)); s1=zeros(size(t));s2=zeros(size(t));
s3=zeros(size(t));s=zeros(size(t));
for p=1:length(t)
for l=1:100
s1(p)=s1(p)+cos(2*pi*l*t(p))*exp(-2*(l*pi)^2*a);
end
end
for p=1:length(t),
for m=1:100,
for n=1:100, 
s2(p)=s2(p)+(1/(m*n))*(cos(2*pi*(m-n)*t(p))*exp(-2*(m-n)^2*c*a)-cos(2*pi*(m+n)*t(p))*exp(-2*(m+n)^2*c*a));
end
end
end
for p=1:length(t),
for m1=1:100
for n1=1:100
s3(p)=s3(p)+(1/(m1*n1))*exp(-(m1^2+n1^2)*c*a)*sin(2*pi*n1*t(p))*sin(2*pi*m1*t(p));
end
end
end
for p=1:length(t),
s(p)=(s1(p)*4+1)*a+(1/c)*(.5*s2(p)-s3(p));
end
plot(t,s)