%espde2.m
ua=0; ub=0;
f=@(x)[-0.1*(abs(x)<0.5)-0.3*(abs(x-5)<0.5)];
xa=-2;xb=10;

 NN=[10,20,40,80,160];
for N=NN
    h=(xb-xa)/(N-1);
x=linspace(xa,xb,N)';
e=ones(N,1);
A=spdiags([-e,2*e,-e],(-1:1),N,N);
A(1,:)=zeros(1,N); A(1,1)=1;
A(N,:)=zeros(1,N); A(N,N)=1;
b=h^2*f(x); b(1)=ua; b(N)=ub;

u=A\b;
 
figure(1)
 plot(x,u,'r')
 legend('numerical')
 xlabel('x')
 axis equal
 axis([-2,10,-1,1])
 end

print(1,'-deps2c','sol_pde2')
