DSolve and NDSolve plots are different

Using FEM will certainly give you the same solution as that of DSolve,

sol2 = NDSolve[{eqn, DirichletCondition[u[x] == 0, x == 100], 
   DirichletCondition[u[x] == 1, x == 0]}, u, {x, 0, 100}, 
  Method -> {"FiniteElement", MeshOptions -> MaxCellMeasure -> 0.001}]

Plot[u[x] /. sol2, {x, 0, 100}, PlotRange -> Full]

My understanding is that when using DirichletCondition, NDSolve will automatically use FEM. So it was using FEM already. The issue seems to be in the grid size used, by default in FEM. since you did not specify one, it used large one and that gave inaccurate result.

k = 1/4;
diffCo = 0.01;
bc = {DirichletCondition[u[x] == 1, x == 0], 
   DirichletCondition[u[x] == 0, x == 100]};
eqn = diffCo*u''[x] + k*u[x] == 0;
solAnalytic = u[x] /. First@DSolve[{eqn, bc}, u[x], {x, 0, 100}];
solNumeric = NDSolve[{eqn == 0, bc},u,{x, 0, 100},Method ->{"FiniteElement"}];

Plot[u[x] /. solNumeric, {x, 0, 100}, PlotRange -> Full]

Mathematica graphics

The smaller the diffCo, the smaller the maxCellMeasure is needed. This manipulate shows this more clearly. When diffCo is large, the default will work OK:

k = 1/4;
diffCo = 1;
bc = {DirichletCondition[u[x] == 1, x == 0], 
   DirichletCondition[u[x] == 0, x == 100]};
eqn = diffCo*u''[x] + k*u[x] == 0;
solAnalytic = u[x] /. First@DSolve[{eqn, bc}, u[x], {x, 0, 100}];
solNumeric = 
  NDSolve[{eqn == 0, bc}, u, {x, 0, 100}, Method -> {"FiniteElement"}];
Plot[u[x] /. solNumeric, {x, 0, 100}, PlotRange -> Full]

Mathematica graphics

Which now matches the DSolve solution.

enter image description here

Added

btw, your ode is of form $\delta u''(x) + k u(x)=0$ where $\delta$ is very small. This calls for boundary layer theory and asymptotic expansion, which is used in these cases, and it gives good result as long as one is close to the point of expansion.

k = 1/4;
diffCo = 1/10^6;
bc = {u[0] == 1, u[100] == 0};
eqn = diffCo*u''[x] + k*u[x] == 0;
solAnalytic = u[x] /. First@DSolve[{eqn, bc}, u[x], x];
Plot[Evaluate@solAnalytic, {x, 0, .05}, PlotRange -> Full, 
 ImageSize -> 300]

Mathematica graphics

(*use 50 terms in expansion around 0 *)
sol = AsymptoticDSolveValue[{eqn, bc}, u[x], {x, 0, 50}];
Quiet[Plot[Evaluate@sol, {x, 0, .05}, PlotRange -> Full, ImageSize -> 300]]

Mathematica graphics

Added Manipulate code

Manipulate[
 Module[{k, u, x, bc, eqn, solAnalytic, solNumeric},
  k = 1/4;
  bc = {DirichletCondition[u[x] == 1, x == 0], DirichletCondition[u[x] == 0, x == 100]};
  eqn = diffCo*u''[x] + k*u[x] == 0;
  solAnalytic = u[x] /. First@DSolve[{eqn, bc}, u[x], {x, 0, 100}];
  solNumeric = NDSolve[{eqn == 0, bc}, u, {x, 0, 100}, 
    Method -> {"FiniteElement", MeshOptions -> MaxCellMeasure -> maxCellMeasure}];

  Grid[{
    {Plot[solAnalytic, {x, 0, 100}, PlotRange -> Full, ImageSize -> 300]},
    {Plot[Evaluate[u[x] /. solNumeric], {x, 0, 100}, 
      PlotRange -> Full, ImageSize -> 300]}
    }, Frame -> All, Spacings -> {1, 1}, Alignment -> Center]
  ]
 ,

 {{diffCo, 0.001, "diffCo"}, 0.001, 1, 0.01, Appearance -> "Labeled"},

 {{maxCellMeasure, 1, "MaxCellMeasure"}, 0.001, 1, 0.001, Appearance -> "Labeled"},
 TrackedSymbols :> {diffCo, maxCellMeasure}
 ]