# Unstable boundary values using NDSolve?

I have an equation:

w1[u_]=-(1/2) u^2 ea e0 Sin[2 a[z]] + K (a''[z])


with numerical values,

K = 6.2*10^-12; ea = 10; e0 = 8.85*10^-12; L = 1;

nd=NDSolve[{w1[1]==0,a[0]==0,a[L]==0},a[z],{z,0,L}];

gr1 = Plot[a[z] /. nd[[1]], {z, 0, L}, PlotRange -> Automatic]


The NDSolve should give the following output for different u, for u equal to zero a is 90 and by increasing u, a decreases till 0. Boundary conditions: @ z=0 and L a is 0 (very known problem),

However, the NDSolve is unstable and not giving the exact solution for the above program. Any suggestions or recommendations.

• what is ele[1]? And where is a[u] used? Jan 10, 2023 at 13:06
• @Nasser made a mistake editing the question.
– a019
Jan 10, 2023 at 13:08
• Well, I doubt that NDSolve is making an error, This tells me you have something wrong with your ode or with some of the numerical values you are using. Are you sure you are using same ode and parameter values used to generate the plot you are showing? Jan 10, 2023 at 13:23
• – a019
Jan 10, 2023 at 13:29
• I had quick look at the paper. Notice it says for fig 4, it is using the normalized z* which is z/h (do not know what h is. Planck's constant?) and your U is E*h. I think the problem is scaling/normalization issue as the ode look OK. !Mathematica graphics May be if you use the normalized values, the plot will show OK since h=6.62607015 × 10-34 which will make big difference in scaling. Jan 10, 2023 at 13:40

This problem can be solved with using the Euler wavelets collocation method described in our paper, we have

w1[u_] = -u^2 Sin[2 a[z]]/kk/2 + (a''[z]); ea = 10; e0 =
8.85*10^-12; kk = 6.2*10^-12/ea/e0; L = 1;
UE[m_, t_] := EulerE[m, t];
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) UE[m, 2^k t - 2 n + 1], (n - 1)/2^(k - 1) <= t <
n/2^(k - 1)}, {0, True}}];
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 7; With[{k = k0, M = M0},
nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; zcol =
xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y; A = Array[as, {nn}];
X2[y_] := A . Psi[y]; X1[y_] = A . int1[y] + c1;
X0[y_] = A . int2[y] + c1 y + c2; eq =
Table[-(1/2) u^2 Sin[2 X0[y]]/kk + X2[y] == 0, {y,
xcol}]; bc = {X0[0] == Pi/2, X0[1] == Pi/2};

var = Join[A, {c1, c2}]; Do[
soln[j] =
FindRoot[Join[eq, bc] /. u -> j,
Table[{var[[i]], 1/10}, {i, Length[var]}]];, {j, 0, 7}]
Plot[Evaluate[Table[2 90 /Pi X0[y] /. soln[j], {j, 0, 7}]], {y, 0, 1},
Frame -> True, PlotLegends -> Table[j, {j, 0, 7}]]



• wonderful stuff!
– bmf
Jan 10, 2023 at 15:52
• @AlexTrounev Can the same method be extended to similar solution but more complex problem such as for two coupled equations and four boundary conditions? The end solution is similar but starting from 2 and stabilizing at 90. I am trying NDSolve , shooting method but after certain range the NDSolve with shooting method can not find the solution. eq1: w1[u_]= L/K3(1/4 ((𝐾2−𝐾3 ) Sin[4 𝜃[z]] 𝜙'[z]^2+2 Sin[2 𝜃[z]] (e0 ea 𝑢'[z]^2+(𝐾3−𝐾1 ) 𝜃'[z]^2+𝐾2 𝜙'[z]^2 ))+(𝐾1 Cos[𝜃[z]]^2+𝐾3 Sin[𝜃[z]] ) (d𝜃'[z])/d𝑧)=0.
– a019
Jan 26, 2023 at 15:17
• @Alex Trounev eq2. w2= L/K3(Cos[𝜃[z]](−2(𝐾2+(𝐾2−𝐾3 ) Cos[2 𝜃[z]] Sin[𝜃[z]] 𝜃'[z] 𝜙'[z]+Cos[𝜃[z]](𝐾2 Cos[𝜃[z]]^2+𝐾3 Sin[𝜃[z]]^2) (d𝜙'[z])/d𝑧))=0. The parameters are now, K1=13.4*10^-12, K2=6*10^-12, K3=18.5*10^-12, ea=11.35, e0=8.85*10^-12, L=1. Boundary conditions, 𝜃[0]== 𝜃[L]==2*(Pi/180), 𝜙[0]==0, 𝜙[L]==Pi/2. And u goes again from 0,to 7. The end solution is similar but starting from 2 and stabilizing at 90. Kind of mirror image of the above. Instead of a[z], now there is 𝜃[z] and 𝜙[z] and {z,0,L}. Extension of the same problem.
– a019
Jan 26, 2023 at 15:20
• @MuhammadAli Yes, it can be extended to solve more complicated problem. Jan 27, 2023 at 4:26
• @AlexTrounev I am not getting the same output, the one you showed above. Is there anything to do with the version of Mathematica I am using (11.0). I am adding a link to the imgur.com/a/6ijQjDU I tried with different things and plotrange but still not.
– a019
Jan 27, 2023 at 8:41