g0 = 1; % Boundary conditions.
gL = 0;

L = 1; % Size of problem domain.
N = 30000; % Number of intervals.
h = L/(N+1); % dx
i = 0:N+1; % Enumerated grid points.
x = (i.*h)'; % x-positions. Column vector.
x_inner = x(2:end-1); % End points are excluded from computations since they are fixed.

f = @(x) x.*exp(x); % r.h.s function.

% Set up matricies
ones_vec = ones(N,1);
A = spdiags([ones_vec, (h^2*cos(x_inner).^2-2).*ones_vec, ones_vec], [-1 0 1],N,N); % See lab2 for more details.
A(1,1) = (h^2*cos(x_inner(1)).^2-1);
full(A);
b = h.^2.*f(x_inner); % b depends on x-position.
% Include the constant endpoint values in boundaries of "inner" grid. Remember, they are
% part of the l.h.s. expression that we have moved over to vector b.
b(1) = b(1) + h;
b(end) = b(end) - gL;
u = A\b; % Solution
u = [u(1)-h; u; gL]; % Add end-points.