clear

%%
% BVP:
%
% -(x^2 * u')' + 4*u = sin(pi*x), 0 <= x <= 1, u(0) = u(1) = 0.
%
% The coefficient functions need to accept a vector input and return a
% vector of the same size. 
c = @(x) x.^2;
s = @(x) 4*ones(size(x));
f = @(x) sin(pi*x);

%%
np = 20;
a = 0;
b = 1;
[x, u] = fem(c, s, f, a, b, np);
[x2, u2] = fem(c, s, f, a, b, 2*np);

plot(x, u, '.-')
hold on
plot(x2, u2, '.-')
xlabel('x')
ylabel('u')
title('Solution by finite elements')
grid on
legend('n points', '2*n points', 'location', 'best')