function u = solve_poisson_jacobi(f, u_exact)

% 1D Poisson equation with Dirichlet boundary conditions

N = length(f) - 1;      % N is a power of 2
h = 1/N;                % true grid: N+1 points
u = zeros(N+1,1);
I = 2:N;
x = (0:h:1)';

for k = 1:50000
    
    u(I) = 1/2 * (u(I+1) + u(I-1) - h^2 * f(I));
    
    if mod(k,10)==0,        % display instructions
            
        figure(1); clf;
        subplot(2,1,1); plot(x,u,'b-'); axis([0 1 -0.1 0.1]);
        err = u - u_exact;
        subplot(2,1,2); plot(x,err,'r-'); axis([0 1 -0.18 0]);
        
    end
    
end
