 x = chebfun('x',[0 1]);
 nmax = 40;
 f = chebfun('1',[0 1]);
 u = (x-x.^2)/2;            % exact solution to -u''(x) = f, u(0)=u(1)=0

 fsum = chebfun('0',[0 1]);
 usum = chebfun('0',[0 1]);

 figure(1), clf
 for n=1:nmax
    lambda_n = (n*pi)^2;
    psi_n = sqrt(2)*sin(n*pi*x);
    cn_cf = psi_n'*f;                     % computed in Chebfun
    cn    = sqrt(2)*(1-(-1)^n)/(n*pi);    % exact integral
    fprintf('step n=%2d; Chebfun computed coefficient: %20.17f\n', n, cn_cf);
    fprintf('step n=%2d; exactly computed coefficient: %20.17f\n\n', n, cn);
    fsum = fsum + cn*psi_n;
    usum = usum + (cn/lambda_n)*psi_n; 
    clf
    subplot(2,1,1)
    plot(f,'b-','linewidth',3);
    hold on
    plot(fsum,'r-','linewidth',3);
    title(sprintf('approximation to $f$ from %d terms in series',n),'fontsize',18,'interpreter','latex') 
    xlabel('$x$','fontsize',18,'interpreter','latex')
    ylabel('$f(x)$','fontsize',18,'interpreter','latex')
    subplot(2,1,2)
    plot(u,'b-','linewidth',3);
    hold on
    plot(usum,'r-','linewidth',3);
    xlabel('$x$','fontsize',18,'interpreter','latex')
    ylabel('$u(x)$','fontsize',18,'interpreter','latex')
    title(sprintf('approximation to $u$ from %d terms in series',n),'fontsize',18,'interpreter','latex') 
    pause
 end
