x_min= 0;
x_max = 1;
N = 10;
L = x_max-x_min;
h = L/N;
x = linspace(x_min,x_max,N+1)'
uexact= @(x) 4*(x-x.^2);
f = @(x) 8*15*ones(1,N-1); % Source term
u_exact = uexact(x);
ua = uexact(x_min);
ub = uexact(x_max);
u = zeros(N-1,1);
A = zeros(N-1,N-1);
b = f(x(2:N)); %
b(1) = b(1) + ua/h^2;
for i=2:N-1
b(i) = b(i) + ub/h^2;
end
dA = diag( 2*ones(1,N-1) );
dAp1 = diag( -1*ones(1,N-2), 1 );
dAm1 = diag( -1*ones(1,N-2), -1 );
A = (dA + dAp1 + dAm1);
A = 15*A/h^2;
u = A\b.'; % solving the linear system
g = [ua; u; ub];
plot(x,u_exact,'b',x,g,'ro-');
legend('exact','numerical');
