function A=slseries(I,u,v,f)
%
%   A = slseries(I,u,v,'f') is a program which enables the user to 
%   examine the expansion of the function f(x) described in f.m in terms 
%   of the eigenfunctions of a Sturm-Liouville problem.  The inputs are 
%   the endpoints of the interval I = [a,b], the matrix u which has the 
%   eigenfunctions evaluated at equidistributed points in the interval, 
%   the matrix v which contains the "dual" functions (which differ from 
%   u only if there is a nontrivial weight function), and the name of 
%   the m-file describing the function.  
%   
%   After computing the vector of coefficients the program presents the 
%   user with the following menu:
%
%   Enter:
%      n   to see the sum of n terms.
%      0   to print the current graphics window.
%     -1   to exit.
%
%   The ouput is a vector containing the coefficients in the
%   eigenfunction expansion.
%

if (nargin ~= 4)
	error('There must be precisely four input arguments.')
end
a = I(1);
b = I(2);
[p,q]=size(u);
disp(['Evaluating your function ...']);
x=(a+[0:p-1]*(b-a)/(p-1))';
fuhjk = zeros(1,p);
for j=1:p
    fuhjk(j) = feval(f,x(j));
end
plot(x,fuhjk);shg
title(['Your function y = ',formula(f),'.']);
A=fuhjk*v;
disp('Enter:');
disp('   n   to see the sum of n terms.');
disp('   0   to print the current graphics window.');
disp('  -1   to exit.');
k=input('');
while k > -1
    if round(k)~=k
	disp(['You must enter an integer between -1 and ',num2str(q),'.']);
    elseif k==0
	print;
    elseif k>q
	disp(['You must enter an integer between -2 and ',num2str(q),'.']);
    else    
	titlest=['The approximation of order ',num2str(k),' to y = ',formula(f),'.'];
	W=u(:,1:k)*A(1:k)';
	plot(x,W,x,fuhjk);shg

	title(titlest);
	E = abs(W(:)-fuhjk(:));
	maxerr = max(E);
	E(1) = (E(1)+E(p))/2;
	l2err = norm(E(1:p-1))*sqrt((b-a)/(p-1));
	xstr=['The maximum error is ', num2str(maxerr), '.    The L2 error is ',num2str(l2err),'.'];
	xlabel(xstr);
    end
    disp('Enter:');
    disp('   n   to see the sum of n terms.');
    disp('   0   to print the current graphics window.');
    disp('  -1   to exit.');
    k=input('');
end