% Parameters
N = 501;             % number of grid points
x = linspace(0,1,N); % grid
dx = x(2)-x(1);      % grid spacing
p_values = [-10, -3, -1, 0, 1, 3, 10];

% Analytical solution function
u_analytical = @(x,p) (abs(p) < 1e-12).*x + (abs(p) >= 1e-12).*((exp(p.*x)-1)./(exp(p)-1));

figure;
hold on;

for k = 1:length(p_values)
    p = p_values(k);
    
    e = ones(N-2,1);
    % Second derivative: -u'' 
    D2 = spdiags([e -2*e e], -1:1, N-2, N-2)/dx^2;
    % First derivative: p*u'
    D1 = spdiags([-e*0.5 e*0.5], [-1 1], N-2, N-2)/dx; 
    L = -D2 + p*D1;
    
    % Right-hand side
    b = zeros(N-2,1);
    
    % Incorporate boundary conditions
    b(1) = b(1) + 0/dx^2;      % u(0)=0
    b(end) = b(end) + 1/dx^2;  % u(1)=1
    
    % Solve linear system
    u_inner = L\b;
    u_num = [0; u_inner; 1]; % add boundary values
    
    % Plot numerical and analytical solution
    plot(x, u_num, '--', 'LineWidth', 2, 'DisplayName', sprintf('num, p=%g',p));
    plot(x, u_analytical(x,p), '-', 'LineWidth', 2, 'DisplayName', sprintf('ana, p=%g',p));
end

xlabel('x');
ylabel('u(x)');
title('Numerical vs Analytical solution of -u'''' + p u'' = 0');
legend show;
grid on;
hold off;