%
%  code24.m -- sturm-liouville greens function (djm - 05 nov 01)
%

N = 20;  tt = pi/2;  TT = 0:pi/(2*N):2*pi;  

%  evaluate exact solution

figure(1);  clf;  subplot(2,1,1)
green = -(TT<tt).*(TT.*(1-tt/2/pi)) -(TT>=tt).*(tt.*(1-TT/2/pi));
plot(TT,green,'g');  hold on;  plot(TT,green,'k.')

%  sum series

errors = zeros(1,N);
coeffs = zeros(1,N);
series = zeros(size(TT));

for j = 1:N
  coeffs(j) = sin(j/2*tt)*4/(j^2)/pi;
  series = series - coeffs(j)*sin(j/2*TT);

  errors(j) = max(abs(green-series)); 
end

%  make plots

plot(TT,series,'b');

title(['\bf Sturm-Liouville Green''s Function (' num2str(N) '-terms)'])  
xlabel('\bf t-axis')
ylabel('\bf G(T;t)=green (what else?) ; g_N(T;t)=blue')
text( pi,-1,'\bf play with the N value')

subplot(2,1,2)
plot(TT,green-series,'b');  hold on;  plot(TT,green-series,'k.');

title(['\bf Error (' num2str(N) '-terms)'])  
xlabel('\bf T-axis')
ylabel('\bf G(T;t)-g_N(T;t)')

figure(2);clf
loglog(1:N,errors,'bx');  hold on
loglog(1:N,abs(coeffs),'ro')

title(['\bf Truncation Error'])  
xlabel('\bf N-axis')
ylabel('\bf max|G(T,t)-g_N(T;t)|=blue ; |coeffs(N)|=red')

plot([3 30],3*[10^(-1) 10^(-2)],'k--')
plot([3 30],5*[10^(-2) 10^(-4)],'k--')

axis([1 N+10 1e-4 1e-0])
text( 5,3e-1,'\bf slope = -1')
text(10,1.e-3,'\bf slope = -2')