NDSolve unable to find initial conditions that satisfy the residual function within specified tolerances

Question

I'm trying to solve a set of differential-algebraic equations with boundary conditions:

$A\frac{d^2\psi}{dx^2}=[c_{F}-2c_{D}(x)+c_{C}(x)]F$

$c_{D}(x)=0.1\exp(-0.2253 - 76.5596\psi(x) - 1.1149)$

$c_{C}(x)=0.2\exp(-1.3305 + 38.2798\psi(x) - 0.3459)$

$\psi\,'(0)=0,\psi(\frac{1}{2})=0$

With the following code, shooting method is used to find the Dirichlet BC at $x=0$ to use NDsolve to solve the DAE

Clear["Global`*"]

A = 5.31251*10^-10;
F = 96485.33252;
cF = 1000;

(*Obtain Dirichlet condition of psi at x=0 by shooting method*)
fun[y_?NumericQ] := Module[{eqns, sol, ics, x, psi, cD, cC},
   eqns = {A*D[psi[x], x, x] == (cF - 2*cD[x] + cC[x])*F,
     cD[x] == 0.1*Exp[-0.2253 - 76.5596*psi[x] - 1.1149],
     cC[x] == 0.2*Exp[-1.3305 + 38.2798*psi[x] - 0.3459]};
   ics = {psi[0] == y, psi'[x] == 0 /. x -> 0};
   sol = 
    NDSolve[{eqns, ics}, {psi[x], cD[x], cC[x]}, {x, 0, 1/2}][[1]];
   psi[x] /. sol /. x -> 1/2];
ic = FindRoot[fun[x], {x, 0}]

dae = {A*D[psi[x], x, x] == (cF - 2*cD[x] + cC[x])*F,
   cD[x] == 0.1*Exp[-0.2253 - 76.5596*psi[x] - 1.1149],
   cC[x] == 0.2*Exp[-1.3305 + 38.2798*psi[x] - 0.3459]};
bcs = {psi[0] == x /. ic[[1]], D[psi[x], x] == 0 /. x -> 0};

(*Solve the differential equation using NDSolve*)
sol = NDSolve[{dae, bcs}, {psi, cD, cC}, {x, 0, 1/2}];

The NDSolve returns the following error:

Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions.

Even if the system of DAE is converted to an ODE to eliminate the need for shooting method, applying NDSolve gives the error:

Clear["Global`*"]

A = 5.31251*10^-10;
F = 96485.33252;
cF = 1000;

ode = {A*
     D[psi[x], x, 
      x] == (cF - 2*0.1*Exp[-0.2253 - 76.5596*psi[x] - 1.1149] + 
       0.2*Exp[-1.3305 + 38.2798*psi[x] - 0.3459])*F};
bcs = {psi[1/2] == 0, D[psi[x], x] == 0 /. x -> 0};

(*Solve the differential equation using NDSolve*)
sol = NDSolve[{ode, bcs}, {psi}, {x, 0, 1/2}];

NDSolve::ndsz: At x == 2.012990504102674`*^-9, step size is effectively zero; singularity or stiff system suspected.

Are there any methods that I can use to get around this? I'm new to numerical method, any help is greatly appreciated, thanks a lot.

Answer 1

Eliminating the unnecessary variables converts the first system from a DAE to an ODE, which gets rid of the problem complained of in the OP:

fun[y_?NumericQ] := 
  Module[{eqns, sol, ics, x, psi, cD, cC}, 
   eqns = {A*
       D[psi[x], x, 
        x] == (cF - 
         2*(0.1*Exp[-0.2253 - 76.5596*psi[x] - 1.1149]) + (0.2*
           Exp[-1.3305 + 38.2798*psi[x] - 0.3459]))*F};
   ics = {psi[0] == y, psi'[x] == 0 /. x -> 0};
   NDSolveValue[{eqns, ics}, {psi[x]}, {x, 0, 1/2}]
   ];

However, the second derivative of psi grows so fast, other problems arise, naturally, such as singularities.

The second derivative as a function of psi:

Plot[(cF - 
    2*(0.1*Exp[-0.2253 - 76.5596*psi[x] - 1.1149]) + (0.2*
      Exp[-1.3305 + 38.2798*psi[x] - 0.3459]))*F/A, {psi[x], 0, 1}]