Matched Asymptotic Expansion Compared with the Numerical Solution

The following equation is "too difficult" for Mathematica to find an analytical solution:

ϵ y'' + 2 y' + (e^y)=0
with y(0)=0 and y(1)=0.

When something is too difficult, Mathematica returns the command that was typed (after saying "running" in the top bar of the window for about 20 seconds).

So we solve the equation by a discrete numerical method, and by a method of matched asymptotic expansions.  The ODE is a two-point boundary-value problem, which we solve using a shooting scheme.

<<Graphics`Colors`

Clear[eps,y,x]

DSolve[{eps*y''[x]+2y'[x]+Exp[y[x]]==0,y[0]==0,y[1]==0},y[x],x]

So we take out pencil and paper and find the following solution using matched asymptotic expansions:

ym[t_,eps_]:=Log[1/2]*Exp[-2 t/eps]-
Log[1/2+t/2]

Plot[ym[t,.01],{t,0,1}]

Note the following property of our approximate solution.  This information will be needed when we attempt to find a numerical solution.

D[ym[t,ep],t]/.t->0//N

To check our solution, we find a numerical solution for a particular value of ϵ.  Mathematica solves ODE's by marching from initial conditions specified at a starting point.  So we need to iterate to satisfy the condition y[1]=1/2 using what is called a "shooting method": we use y'[0]=(s/ϵ) where we expect s will be a number near 1.

eps=.01;

nshoot[eps_,s_]:=NDSolve[{eps*y''[t]+2*y'[t]+Exp[y[t]]==0,y[0]==0,y'[0]==s/eps},{y},{t,0,1}]

tryit=nshoot[eps,1.]

Let's try to use that InterpolatingFunction thing.  The replacement rule for y is in a nested list.  Let's evaluate it at t=0.3

y[.3]/.tryit[[1,1]]

Plot[y[t]/.tryit[[1,1]],{t,0,1}]

As seen above, with s=1, our shooting scheme missed the boundary condition y[1]=(1/2).  So we need to find the correct value of s.

yAt1[eps_,s_]:=y[1]/.nshoot[eps,s][[1,1]]

yAt1[eps,1.0]

yAt1[eps,2]

Plot[yAt1[eps,s],{s,1,2}]

theroot=FindRoot[yAt1[eps,s]==0,{s,{1.,2.}}]

yAt1[eps,s/.theroot]

So we see that we have a numerical solution that satisfies both boundary conditions. Let's plot it out and compare it with our approximate solution.

numsol=nshoot[eps,s/.theroot]

y[.3]/.numsol[[1,1]]

ym[.3,eps]

Next we plot both the numerical and the perturbation solution together.  They very nearly overlap:

Plot[{y[t]/.numsol[[1,1]],ym[t,eps]},{t,0,1},PlotRange->{0,1},PlotStyle->{Green,Red}]

ydiff[t_]:=(y[t]/.numsol[[1,1]])-ym[t,eps]

ydiff[.3]

Plot[ydiff[t],{t,0,1}]