How to solve the continuity equation for dark and light scenario by symbolic function with boundary condition

Now I met a problem on my research. I met with two ODEs with identical form in left side but one is homogenous and the other is nonhomogeneous.

Here is my code:

syms V_T mu_p F r_m           % constant of ODE1
syms N_v phi_3 phi_4 d_j      % constants for boundary condition
syms p(x) G(x)                        % p means hole density G(x) means generation rate
eqn = -V_T*mu_p*diff(p,x,2) + F*mu_p*diff(p,x) + r_m*p == 0;   %MARKED
P_dark = dsolve(eqn) %% no ; so the form is readable
eqn = -V_T*mu_p*diff(p,x,2) + F*mu_p*diff(p,x) + r_m*p == G(x);   %MARKED
P_ph = dsolve(eqn)                                                          
P_ph_int = int(P_ph,[350*10^-6 800*10^-6])                     % need to integrate the P_ph in the interval of wavelength
P = P_dark+P_ph;
cond = [p(0) == N_v * exp(-phi_3/V_T), p(d_j) == N_v * exp(-phi_4/V_T)]; 
%This is the problem. I know how to handle it if the boundary is for P_ph or P_dark only. But right now, the boundary condition is used for P, the addition of the P_dark and P_ph. I can get the solution with undetermined constants but I can not use the boundary condition to find such constants.
P = dsolve(P,cond)

the bug info is

Error using mupadengine/feval (line 195)

No differential equations found. Specify differential equations by using symbolic functions.

Error in dsolve>mupadDsolve (line 340)

T = feval(symengine,'symobj::dsolve',sys,x,options);

Error in dsolve (line 194)

sol = mupadDsolve(args, options);

Answer 1

[Walter Roberson]

syms V_T mu_p F r_m           % constant of ODE1
syms N_v phi_3 phi_4 d_j      % constants for boundary condition
syms p(x) G(x)                        % p means hole density G(x) means generation rate
eqn = -V_T*mu_p*diff(p,x,2) + F*mu_p*diff(p,x) + r_m*p == 0;  % equation in dark 
P_dark(x) = dsolve(eqn) 

P_dark(x) = 

eqn = -V_T*mu_p*diff(p,x,2) + F*mu_p*diff(p,x) + r_m*p -G(x) == 0;  % equation in light (your suggestion)  %MARKED
P_ph(x) = dsolve(eqn)

P_ph(x) = 

new_vars = setdiff(symvar(P_ph(x)), symvar(eqn))

new_vars = 

P_total(x) = P_dark+P_ph  % Solutions which boundary condition applied

P_total(x) = 

eqn1 = P_total(0) == N_v * exp(-phi_3/V_T)  %% boudary condition of P at x=0

eqn1 = 

eqn2 = P_total(d_j) == N_v * exp(-phi_4/V_T) %% boudary condition of P at x=d_j

eqn2 = 

eqns = [eqn1;eqn2]

eqns = 

symvar(eqns)

ans = 

S = solve(eqns, new_vars)

S = struct with fields:
    C1: [1×1 sym]
    C2: [1×1 sym]

temp = struct2cell(S);
new_vars(:) == vertcat(temp{:})

ans = 

[Jiaxing Liu]

Hi Walter,

Thanks for your answer! You solved my problem.

I ran your code. But I find some differences in the resutls.

You got the same undetermined constant C_1 C_2 in both dark and light scenario.However, I get different constants for these two conditions. I think because the two equation have different layout. There should be four constant such as C_3 C_4 C_5 C_6.

So My question is:

1. I wonder if you have some idea on the reason of this problem?

Here is the results I ran your code in R2019a

P_dark(x) =

P_ph(x) =

Answer 2

[Torsten]

The equation for P reads

-V_T*mu_p*P''+ F*mu_p*P'+ r_m*P - G(x) = 0.

Now you can apply the boundary conditions.

[Jiaxing Liu]

Hi Torsten,

Thanks for your advice.

I try your answer but it seems that it don't change my situation.

I think the problem is how to solve the equation of P (which is the addition of the solution of P_dark and P_ph_int). So the problem should be in the last row of the codes.

So I wonder if you could help me on finding ways to solve the P using solve function with boundary conditions?

Jiaxing Liu

[Torsten]

In your code from above, P = P_dark + P_ph, not P = P_dark + P_ph_int.

What do you get as error message if you try to solve "my" differential equation for P with your boundary conditions ?

[Jiaxing Liu]

Hi Torsten,

Thanks!

I used your differential equation and I change the way to not do the integration. Because the integration can be done before the equation. It will make the life easier.

I think the problem has been solved with new questions comes.

syms V_T mu_p F r_m           % constant of ODE1
syms N_v phi_3 phi_4 d_j      % constants for boundary condition
syms p(x) G(x)                        % p means hole density G(x) means generation rate
eqn = -V_T*mu_p*diff(p,x,2) + F*mu_p*diff(p,x) + r_m*p == 0;  % equation in dark 
P_dark(x) = dsolve(eqn) 
eqn = -V_T*mu_p*diff(p,x,2) + F*mu_p*diff(p,x) + r_m*p -G(x) == 0;  % equation in light (your suggestion)  %MARKED
P_ph(x) = dsolve(eqn)                                                          
P_total(x) = P_dark+P_ph  % Solutions which boundary condition applied
eqn1 = P_total(0) == N_v * exp(-phi_3/V_T)  %% boudary condition of P at x=0
eqn2 = P_total(d_j) == N_v * exp(-phi_4/V_T) %% boudary condition of P at x=d_j
eqns = [eqn1,eqn2]
S = solve(eqns)

After did in this way, the code can run without error. But I can not get the undetermined constant I want. So the new problem comes that "how to use solve function to get the undertermined constant?"

Jason