function mit18336_poisson1d_error
%MIT18336_POISSON1D_ERROR
%    Sets up and solves sequence of 1d Poisson problems on
%    domain [0 1], with homogeneous Dirichlet b.c.
%    Plots the convergence of the error in log-log scale.

% 02/2009 by Benjamin Seibold
%            http://www-math.mit.edu/~seibold/
% Feel free to modify for teaching and learning.

n = ceil(10.^linspace(1,4,20));                % sequence of cell numbers
h = 1./(n+1);                            % corresponding grid resolutions
e = n*0;
clf
for i = 1:length(n)
    x = (h(i):h(i):1-h(i))';       % regular grid without boundary points
                                                      % 1D Poisson matrix
    K1D = spdiags(ones(n(i),1)*[-1 2 -1],-1:1,n(i),n(i))/h(i).^2;
    f = f_corr(x);                                      % right hand side
    u = K1D\f;                                      % solve linear system
    u_diff = u-u_corr(x);                          % calculate difference
    u_diff = u_diff/norm(u_diff,inf);
    subplot(1,2,1), plot(x,u_corr(x),'b-',x,u,'r-',x,u_diff,'k-')
    legend('true solution','numerical approximation','scaled difference')
    title('1d poisson equation')
    subplot(1,2,2)
    e(i) = norm(u-u_corr(x),inf);           % calculate error in max norm
    loglog(10.^[-4 0 -4],10.^[-4 0 -8],'k-',h,e,'b.-')
    title('error convergence')
    drawnow
end

%========================================================================

function y = u_corr(x)                                 % correct solution
c = 9*pi;
y = sin(c*x.^2);

function y = f_corr(x)                                  % right hand side
c = 9*pi;
y = sin(c*x.^2).*(2*c*x).^2-cos(c*x.^2)*(2*c);