%Solution of the model problem BVP with FDM, constant stepsize h
%-u"=f(x), u(0)=u(1)=0, f(x)=sin(pi*x)
a=0;b=1;  %interval bounds
N=63;  %N+1=number of steps
h=(b-a)/(N+1);  %stepsize
xp=a:h:b;  %DG, discretize interval to a grid
x=a+h:h:b-h; %inner points generated
f=h*h*[sin(pi*x)]'; %right hand side generated
A=zeros(N,N); 
for i=1:N-1
    A(i,i)=2;
    A(i,i+1)=-1;
    A(i+1,i)=-1;
end
A(N,N)=2;  %full matrix A generated
u=A\f;   %solution of linear system of equations
up=[0,u',0]; %add the boundary values to the solution vector
plot(xp,up)
title('Solution of -u"=sin(pi*x), u(0)=u(1)=0')
xlabel('x')
ylabel('u')
grid