function [x,f] = findho(n)

% -----------------------------------------
%
% Computes the first excited state of the harmonic 
% oscillator
% 
% Lf = lam*f, here lam = 3 
%
% where L = -d^2/dx^2+x^2 with Dirichlet b/c
% 
% Note: n needs to be odd
%
% Example: 
%
% [x,f] = findho(1001);
% y = (x.*exp(-x.^2/2))';
% a = max(f)/max(y);
% plot(x,f-a*y)
%
% ------------------------------------------

lam = 3;

A = sparse(n+1,n);

x = linspace(-10,10,n);

dx = x(2)-x(1);
dxs = dx^2;

myind = floor(n/2+1/dx);

A(1,myind) = 1;

A(2,1) = 2/dxs+x(1)^2-lam;
A(2,2) = -1/dxs;

A(2,n) = 1/dxs;

for j=3:n
    A(j,j-2) = -1/dxs;
    A(j,j-1) = 2/dxs+x(j-1)^2-lam;
    A(j,j) = -1/dxs;
end

A(n+1,1) = 1/dxs;

A(n+1,n-1) = -1/dxs;
A(n+1,n) = 2/dxs+x(n)^2-lam;

b = zeros(n+1,1);

b(1) = 1;

f = A\b;

end