Solving two second order BVPs and the boundary conditions are related

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));
%-----------------------------------------------------------------------%

[darova]

Maybe solve 2 equations separately?

[Mohammed Elhefnawy]

I tried but how to connect the boundary conditions of the two equations?

[darova]

Solve for region 2 (-b .. 0)

Get solution and solve for region 1 (0 .. a)

[Mohammed Elhefnawy]

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

Answer 1

[darova]

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]

(I don't know if this can help) the analytical solution is (before applying BCs) is

[darova]

How to compute A,B,C,D constants?

[Mohammed Elhefnawy]

I think by solving the four equations formed when substituting by the boundary conditions in

(u1 &u2) equations after solving them

Thanks alot for you time and effort

I really appreciate it

[darova]

No matter what initial conditions i try the results always looks the same (u1=u2)

I found that k,A and B are similar. Do you have original formulas for these?

k=(2*m*E)/(hbar^(2));
V=-(1.6*(10^(-19)))/(4*3.14*8.854*(10^(-12))*abs(x));    % V is very small
A=sqrt((2*m*E)/(hbar^(2)));
B=sqrt(((2*m)/(hbar^(2)))*(E-V));

Try

b = 5e-42;
a = 5e-42;
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*1e-19/(4*3.14*8.854*1e-12*abs(x));
A = 2*m*E/hbar^2;
B = @(x) 2*m/hbar^2*(E-V(x));
myode = @(x,u) [...
    u(2)
    -2i*k*u(2)+(k-A)*u(1)
    u(4)
    -2i*k*u(4)+(k-B(x))*u(3)];

u0 = [0 1 0 1]*1e-2;            % initial conditions

[x,u] = ode45(myode, [-b a], u0);

plot(x,u(:,1),'-b')
hold on
plot(x,u(:,3),'.r')
hold off
legend('u1','u2')

[Mohammed Elhefnawy]

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]

(always). Ok?

[Mohammed Elhefnawy]

Yes,rigth.