function a = wrml(xmax,N);
%WRML    This program computes the Least Squares approximation for f(x), where
%
%             f'(x) - f(x) = 0,   f(0) = 1
%
%        on the interval [0,xmax]. The approximate solution is 
%             
%             g(x) = 1 + a(1) x + a(2) x^2 + ... a(N) x^N.
%
%        a = wrml(xmax,N) returns the vector a with xmax and N as inputs.
%

%        Ellen McGrattan, 8-28-96


C     = eye(N) + diag(2:N,-1);
k     = [2:N+1];
j     = k'*ones(1,N)+ones(N,1)*[1:N];
ixx   = (xmax.^j)./j; 
e     = [1;zeros(N-1,1)];
xk    = xmax.^k./k;
iex   = e*xk;
iee   = xmax*e*e';
A     = C*ixx*C'+iex*C'+C*iex'+iee;
b     = C*xk'+xmax*e;
a     = -A\b;
A,b
