% HW 8, Bonus 3
% The only difference from Problem 5 is that the BC at the right end
% is treated with 1st order accuracy. 

clear all

xspan = [0 1];

% boundary conditions
alfa = 0;
beta = 2;
A1 = 1;
A2 = 2;

h = 0.1/2;
xplot = xspan(1) :h: xspan(2);
x = xplot(2:end-1);
N = length(x);

Q = -2./(1+x).^2;
R = -4./(1+x).^2;

c_vec = ones(1,N-1);
a_vec = c_vec;
b_vec = -(2 - h^2*Q);
b_vec(end) = b_vec(end) - A2/(-A2-h*A1);
r_vec = h^2*R;
r_vec(1) = r_vec(1) - alfa;
r_vec(end) = r_vec(end) + h*beta/(-A2-h*A1);
Y = thomas(a_vec, b_vec, c_vec, r_vec);
Yend = (-A2*Y(end) - h*beta)/(-A2-h*A1);
Y = [alfa Y Yend];

Yexact = 2*xplot./(1+xplot);

figure(81201);
subplot(211)
plot(xplot,Y, 'b', xplot,Yexact, 'r--', 'linewidth',2)
xlabel('x', 'fontsize', 12);	ylabel('y', 'fontsize', 12);
title(['Finite-difference method (h = ',num2str(h),')'],'fontsize',11) 
legend('numerical','exact','Location','Best')
subplot(212)
plot(xplot,Y-Yexact)
xlabel('x', 'fontsize', 12);   ylabel('error', 'fontsize', 12)


