function [E,V] = anharmosc(ndim,kappa,lambda,wieviel)

% Berechnung des Eigenwertspektrums des anharmonischen Oszillators
% ndim Dimension der Matrix
% wieviel < ndim Anzahl der zu berechnenden Eigenl"osungen
% Beispiel [E,V] = anharmosc(100,0.5,1,5);

% W. Schweizer
% Vers. 1.1 06/2000 
% Vers. 1.2 03/2016
%
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016


% Preallokation
Nsqrt=zeros(ndim);
E0 = Nsqrt;
Emat = Nsqrt;
x3 = Nsqrt;
x4 = Nsqrt;

qz = [0:ndim-1].';
N = diag(qz);
Nsqrt = sqrt(N);

E0 = N + 1/2*eye(ndim);

ad = [Nsqrt(:,2:end),zeros(ndim,1)];
a = ad';

x = sqrt(1/2)*(ad + a);

x3 = x^3;
x4 = x3*x;

Emat = E0 + 1/8*sqrt(2)*kappa*x3 + 1/8*lambda*x4;

if wieviel >= ndim
    wieviel = ceil(ndim/2),
end

[V,E]=eigs(Emat,wieviel,'sm');
E = diag(E);
spy(Emat),shg
nnz(Emat)/ndim^2

eigenvisu(V,ndim)

function eigenvisu(V,ndim)

% Berechnung der Eigenfunktionen des harmonischen Oszillators
% sqrt(hbar/m omega) = 1

x = linspace(-3,3,2000)'; % normierte Ortskoordinate

wellfunho = [];
for n=0:ndim-1 
    nfak = 1/sqrt(sqrt(pi)*2^n*factorial(n));
    Hn = polyval(hermitpoly(n),x);
    wellfunho = [wellfunho,nfak*Hn.*exp(-x.^2)];
end

wellfunah = wellfunho*V;
%size(wellfunah)
pos = round(length(x)/2);

% Visualisierung
for n=1:size(V,2)
    figure
    vor=sign((wellfunah(pos,n)-wellfunah(pos+10,n))*(wellfunho(pos,n)-wellfunho(pos+10,n)));
    if vor == 0;
        vor =1;
    end
    plot(x,vor*wellfunah(:,n),x,wellfunho(:,n)),grid on,shg
    title([num2str(n),'. Wellenfunktion'])
    legend('anharmon','harm')
end
