function poissonDriver()
    %Wipe memory
    %clc;
    clf;
    clear;

    %Discretization parameters
    N = 2^4;
    xmin = 0;
    xmax = 1.0;
    L = xmax-xmin;
    dx = L/(N-1);
    %Generate mesh
    X = xmin:dx:xmax;
    Uex = Uexact(X);
        
    %Calculate Poisson matrix
    A = getPoisson(N,dx);
    b = getRHS(X)';
    
    %Rewrite Dirichlet BC's
    A(1,:) = 0;
    A(N,:) = 0;
    A(1,1) = 1;
    A(N,N) = 1;
    b(1) = 0;
    b(N) = 0;
    As = sparse(A);
    %Solve discrete system with different solvers
    tic
    fprintf('Solve with inverse\n');
    uh = inv(A)*b;
    toc
    
    tic
    fprintf('Solve with backslash\n');
    uh = A\b;
    toc
    
    tic
    fprintf('First convert to sparse format, then solve with inverse\n');
    uh = inv(As)*b;
    toc
    
    tic
    fprintf('First convert to sparse format, then solve with backslash\n');
    uh = As\b;
    toc
    
    %Plot result and exact solution
    plot(X,uh,X,Uex,'o')

    %Compute error
    err = sqrt(dx*(Uex-uh')*(Uex-uh')');
    fprintf('dx: %e error: %e\n',dx,err);
    fprintf('Done!\n');
end

function y = Uexact(X)
    y = sin(pi*X);
end

function y = getRHS(X)
    y = -pi*pi*sin(pi*X);
end

function A = getPoisson(N,dx)
    A = zeros(N,N);
    for i = 1:N
        A(i,i) = -2;
    end
    for i = 1:N-1
        A(i,i+1) = 1;
        A(i+1,i) = 1;
    end
    invh2 = 1./(dx*dx);
    A = A*invh2;
end