% espde1.m
xa=0;xb=10;
uex=@(x)sin(x)./exp(x); ua=uex(xa); ub=uex(xb);
f=@(x)2*cos(x)./exp(x);


Err=[];
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;
Uex=uex(x);
 plot(x,u,'r',x,Uex,'b')
 legend('numerical','exact')
Err=[Err;max(abs(u-Uex) )];
if N==40
     print(1,'-deps2c','sol_pde1_40')
end
pause
end

figure(2)
H=(xb-xa)./(NN-1);
loglog(H,Err,H,H.^2)
xlabel('h');
grid on
legend('error','h^2')
   print(2,'-deps2c','err_pde1_40')