function sol=harmoscbvp(a,ty)

% Beispiel: Quantenoszillator im Kasten
%
% Aufruf sol = harmoscbvp(2.5,'asym');
%
% Berechnung der quantenmechanischen L"osung
% a ist die Breite des Potentials
% ty: sym   positive Paritaet
%     asym  negative Paritaet


% W. Schweizer
% Vers. 1.1 06/2008 - 01/2016
%
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016

% bvpinit
switch ty
    case 'sym'
           lambda = 0.5;
           symf=1;
           solinit = bvpinit(linspace(0,a,10)...
                    ,@oscinit,lambda,a,ty);
    case 'asym'
          symf=-1;
          lambda = 1.5;
          solinit = bvpinit(linspace(0,a,10)...
                    ,@oscinit,lambda,a,ty); 
end

sol = bvp4c(@oscode,@(ya,yb,lambda) oscbc(ya,yb,lambda,ty) ,solinit);

% Ausgabe und Visualisierung
fprintf('Der Eigenwert ist %7.3f.\n',...
        sol.parameters)

xint = linspace(0,a);
Sxint = deval(sol,xint);
xtot = [-fliplr(xint(2:end)), xint];
ytot = [symf*fliplr(Sxint(1,2:end)),Sxint(1,:)];

plot(xtot,ytot), hold on
ygrs = get(gca,'YTickLabel');
ygr(1) = str2double(ygrs(1,:));
ygr(2) = str2double(ygrs(end,:));
plot([0,0],ygr,':')
ht=title('Eigenfunktion harmonischer Oszillator im Kasten: Abb. 11.2');
ht.FontSize = 9;
ht.FontWeight = 'normal';
shg
hold off

% ------------------------------------------------------------
function dydx = oscode(x,y,lambda,varargin)
% Schroedingergleichung
dydx = [  y(2)
         -(2*lambda - x.*x)*y(1) ];

     % ------------------------------------------------------------
function res = oscbc(ya,yb,lambda,ty)
% Randbedingungen
switch ty
    case 'sym'
           res=[ya(1)-1
                yb(1)
                ya(2)];
    case 'asym'
           res=[ya(1)
                yb(1)
                ya(2)-1];
end

% ------------------------------------------------------------
function yinit = oscinit(x,a,ty)
% Schaetzer
switch ty
    case 'sym'
           yinit=[exp(-1/2*x*x)
                  -x*exp(-1/2*x*x)];
    case 'asym'
           yinit=[x*exp(-1/2*x*x)
                  -x*x*(exp(-1/2*x*x)-1)];
end
