function [u]=diffin(N)
%
h = 1/(N+1);
coeff=(N+1)^2;
I=speye(N,N);diag = 2*I;subd=-sparse(2:N,1:N-1,1,N,N);
A=coeff*(diag+subd+subd');
b=sparse(N,1);
for i=1:N
  b(i) = f(i*h);
end
u=A\b;
%
%   error
%
err=sparse(N,1);
for i=1:N
   err(i)=u(i)-uexact(i*h);
end
fprintf('err inf=%e \n',norm(err,inf))
fprintf('err 2=%e \n',norm(err,2))
fprintf('err A=%e \n',sqrt(dot(err,A*err)))
end
%
% rhs of -u''(x)=f(x)
%
function f = f(x)
%f = 2;
f = sin(2*pi*x)*4*pi*pi;
end
%
% exact sol of -u''(x)=f(x), u(0)=0=u(1)
%
function uexact = uexact(x)
%uexact = x*(1-x);
uexact = sin(2*pi*x);
end