function [ x_out, y_out ] = jacobi_ivp( xmin, xmax, y0 )
   N = size(y0,1);
   h = (xmax - xmin) / (N-1.0);
   error = 1;
   
   while (error > 1e-8)
       % internal difference stencil
       L_y = (-2*y0 + circshift(y0,1) + circshift(y0,-1))/(h*h) + exp(y0);
       D = ones(N,1)*((-2.0)/(h*h)) - exp(y0);

       % endpoints
       L_y(1) = y0(1);
       L_y(end) = y0(end);

       D(1) = 1;
       D(end) = 1;
       error = norm(L_y);
       
       y0 = y0 - L_y./D;
   end
   
   x_out = linspace(xmin, xmax, N),
   y_out = y0;
end

