function [u, f] = lp_exact1d(kind, x)
% LP_EXACT1D  A manufactured solution of -u'' = f on [0, 1], and its source.
%
% Pick u first, differentiate it by hand, and call the result f.  Then u is by
% construction the exact solution of the boundary value problem
%
%     -u''(x) = f(x)  on (0, 1),   u(0), u(1) given by u itself,
%
% and any difference between a computed solution and u is error, known to full
% precision at every point and at no cost.  The five choices below differ only
% in how hard they are to resolve.
switch kind
    case 'smooth'
        % entire, moderate amplitude: the easy case
        s = sin(3*pi*x);
        u = exp(s);
        f = (3*pi)^2 * (s - cos(3*pi*x).^2) .* u;
    case 'oscil'
        % eight wavelengths across the interval: still entire, but a mesh must
        % resolve the wavelength before any convergence rate appears
        u = sin(8*pi*x) + x;
        f = (8*pi)^2 * sin(8*pi*x);
    case 'peak'
        % a Gaussian of width 1/20 at x = 0.4: entire, but with all of its
        % structure in a tenth of the interval
        a = 200;
        u = exp(-a*(x - 0.4).^2);
        f = (2*a - 4*a^2*(x - 0.4).^2) .* u;
    case 'layer'
        % a boundary layer of width eps at x = 1, as in a singular perturbation
        ep = 0.02;
        u = exp(-(1 - x)/ep);
        f = -u/ep^2;
    case 'lowreg'
        % u is C^2 but not C^3: |x - 1/2|^(5/2) has a square root in u''.
        % Every polynomial-based method feels this.
        c = 0.5;
        u = abs(x - c).^(5/2);
        f = -(15/4) * sqrt(abs(x - c));
    otherwise
        error('lp_exact1d: unknown solution "%s"', kind);
end
end