% pset2_1.m 
% Please replace question marks (???) by appropriate code
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=???; % appropriate value of n so delta functions fall on mesh points
h=1/(n+1); % step size
x=0:h:1-h;x=x'; % mesh points
xdelta=[1/2 2/3]; % contains points at which there is a delta function

%%%%%%%%%%%%%%%%%%%%%%%%%%% Finite differences %%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=sparse([],[],[],n+1,n+1,3*(n+1));
T=spdiags(ones(n+1,3)*diag([-1 2 -1]),[ -1 0 1 ],T);T(1,1)=1;
f=zeros(n+1,1); f(round(xdelta/h)+1)=(???);% right-hand side, ONE MORE UNKNOWN
f(1)=f(1)/2; % fix right-hand side for second order accuracy at boundary


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u=(1/h^2)*T\f;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Error %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following creates the true solution, piece-wise
% The code (a<=x)&(x<b) creates a vector the size of x, with 0 everywhere
% except where BOTH a<=x AND x<b.
ut=(???).*((  0<=x)&(x<1/2))+...
   (???).*((1/2<=x)&(x<2/3))+...
   (???).*((2/3<=x)&(x<1  ));
e=sqrt(h*sum(abs(u-ut).^2));disp(e); % L2 error in interval [0,1]
