%
%  lect3.m -- code for BVP example
%
%      y" - s^2 y = x    on x = [0,1]  
%
%      y(0) = y(1) = 0
%
%  finite difference method
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  set parameters 

N = 64;  h = 1/N;
s =  1;  c = - 2 - (s*h)^2;

%  construct diagonal matrices

T = diag(ones(N-2,1),-1) + c*diag(ones(N-1,1),0) + diag(ones(N-2,1),1);

%  x-axis vector

Xj   = [h:h:1-h]';
X    = [0 ; Xj ; 1]; 

% y-axis vectors (including op count for mtx solve)

flops(0)
Yj   = T\(h^2*Xj);
fl = flops;

disp([' '])
disp(['   operation count = ' num2str(fl)])

%  computed & exact solution

Ybar = [0 ; Yj ; 0];
Yex  = sinh(s*X)/(s^2 * sinh(s)) - X/(s^2);

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
%  solution plot

clf;subplot(2,1,1)

plot(X,Ybar,'.')
hold on
plot(X,Yex,'k')

title('my first numerical BVP')
xlabel('x-axis')
ylabel('y-axis')

text(0.4,-0.02,['discretization, N = ' num2str(N)])

%  error plot

subplot(2,1,2)

e = abs(Ybar-Yex);  emax = max(e);
disp([' '])
disp(['   max error = ' num2str(emax)])
disp([' '])

plot(X,e,'.')

title('errors:  e(x) = |y_{exact} - y_{comp}|')
xlabel('x-axis')
ylabel('e-axis')

text(0.4,emax/3,['max error = ' num2str(emax)])