# Solving two second order BVPs and the boundary conditions are related

2 views (last 30 days)
Mohammed Elhefnawy on 13 Dec 2019
I am trying to sove these two BVP equations (solving for time independent schrodinger equation in 1D)
the bounadry conditions are related and
the BCs are
and the analytical solution is (before applying BCs)
I am getting some errors which are shown in the code below
I would appreciate any help please
function dudx = osc(x,u)
h=6.62606896E-34; %% Planck constant [J.s]
hbar=h/(2*pi);
q=1.602176487E-19; %% electron charge [C]
%m=input('enter the mass ');
m=9.10938188E-31; %% electron mass [kg]
%E=input('enter the energy ');
E=10000;
k=(2*m*E)/(hbar^(2));
V=-(1.6*(10^(-19)))/(4*3.14*8.854*(10^(-12))*abs(x));
A=sqrt((2*m*E)/(hbar^(2)));
B=sqrt(((2*m)/(hbar^(2)))*(E-V));
dudx = [u(2)
-2*1i*k*(u(2))+((k^(2))-(A^(2)))*u(1)
u(4) %% error :index exceeds array bounds
2*1i*k(u(4))+((k^(2))-(B^(2)))*u(3)]; %% error: Array indices must be positive integers or logical values.
end
%-----------------------------------------------------------------------%
function res = bcfcn(ya,yb,yc) % boundary conditions
res = [ya(1)-ya(3)
ya(2)-ya(4) %%error : the function bcfcn should return a column of two vectors
yb(1)-yc(1) %% and I am m=not sure if it is right to write the BCs like this
yb(2)-yc(2)];
end
%-----------------------------------------------------------------------%
function g = guess(x) % initial guess for y and y'
g = [exp(x)+exp(-x)
exp(x)+exp(-x)];
end
%-----------------------------------------------------------------------%
xmesh = linspace(-10,10,21);
solinit = bvpinit(xmesh, @guess);
sol = bvp4c(@osc, @bcfcn, solinit);
plot(x,u(:,1));
%-----------------------------------------------------------------------%

Show 1 older comment
Mohammed Elhefnawy on 14 Dec 2019
I tried but how to connect the boundary conditions of the two equations?
darova on 14 Dec 2019
Solve for region 2 (-b .. 0)
Get solution and solve for region 1 (0 .. a)
Mohammed Elhefnawy on 14 Dec 2019
Thanks alot
but there is something that I still don'y get it
I will solve for the second region but to solve it I have to put boundary conditions
so if I assume for example that (u2(0))=some constant (A) and the same for the other conditions and solve it
how to take the solution and pass it to the other region and how to deal with the BCs their too

darova on 14 Dec 2019
How to set up initial guess? Any ideas?
function main
%-----------------------------------------------------------------------%
b = 10;
a = 10;
h=6.62606896E-34; %% Planck constant [J.s]
hbar=h/(2*pi);
q=1.602176487E-19; %% electron charge [C]
%m=input('enter the mass ');
m=9.10938188E-31; %% electron mass [kg]
%E=input('enter the energy ');
E=10000;
k=(2*m*E)/(hbar^(2));
V = @(x) -(1.6*(10^(-19)))/(4*3.14*8.854*(10^(-12))*abs(x));
A = sqrt((2*m*E)/(hbar^(2)));
B = @(x) sqrt(((2*m)/(hbar^(2)))*(E-V(x)));
myode = @(x,u) [u(2)
-2*1i*k*(u(2))+((k^(2))-(A^(2)))*u(1)
u(4)
-2*1i*k*(u(4))+((k^(2))-(B(x)^(2)))*u(3)];
x0 = fsolve(@F,[1 1 1 1]*exp(-b)); % what initial guess should be set?
[x1,u11] = ode45(myode, [-b a], x0);
plot(x1,u11)
function res = F(x0)
[x,u] = ode45(myode, [-b a], x0);
u1 = u(:,1);
du1 = u(:,2);
u2 = u(:,3);
du2 = u(:,4);
res = [interp1(x,u1,0) - interp1(x,u2,0) % u1(0) == u2(0)
u1(end) - u2(1) % u1(a) == u2(-b)
interp1(x,du1,0) - interp1(x,du2,0) % du1(0)== du2(0)
du1(end) - du2(1)]; % du1(a)== du2(-b)
end
end

Mohammed Elhefnawy on 15 Dec 2019
A is (alpha) and B is(Beta)
all of this is to try solving time independent schrodinger equation
these equations are form the book (Semiconductor physics and devices by donald A.Neeman)
chapter 3 section 1
darova on 15 Dec 2019
(always). Ok?
Mohammed Elhefnawy on 15 Dec 2019
Yes,rigth.