function [Z,f] = gaussZ(N,c,delta)
%function [Z,f] = gaussZ(N,c,delta)

B = 0;k=1;

for Z=0:0.1:c
    
    for n=1:1:N        
        Qn = quadl(@myfun,0,c,1.e-9,0,c,n,delta);
        B = B + Qn*sin(n*pi*Z/c);
    end
    
    f(k) = B;
    k = k+1;
    B = 0;
end
Z=0:0.1:c;
        
        
