Precision in Integrals and NDSolve

improved answer

Here I'll show how to solve your problem which is singular for x==0 and x==1 and illscaled for N1>>1

First the differential equation

s0 = 0.01 + 0.007/1000 Sin[th] // Rationalize;
h = 0.3 // Rationalize;
th = Pi/4;
ode[N1_] := (x (1 - x))/(2 N1) D[T0[x], x, x] +  
             s0  x (1 - x)  (x + h (1 - 2 x )) D[T0[x], x] == -Integrate[
         E^(-  N1 s0 a (a + 2 h (1 - a))), {a, 0, x}]/
       Integrate[E^(- N1 s0 a (a + 2 h (1 - a))), {a, 0, 1}] 

solution of the ode N1=10000

T = With[{e = 10^-10, N1 = 10000}, 
NDSolveValue[{ode[N1],
         T0[e] == 0, T0[1 - e] == 0}, {T0}, {x, e, 1 - e} ][[1]] ] 

Plot[T[x], {x, 0, 1}] 

That solves your problem in principle.

I don't have more time to elaborate the case N1=100000 which needs numerical finetuning concerning WorkingPrecision&Co ...

case N1=100000

Thanks to the answer @Dominic ,which let me to the idea using $MaxExtraPrecision = 100 . Blockinstead of With (don't know why it's necessary) evaluates to

T = Block[{$MaxExtraPrecision = 100, e = 10^-10, N1 = 100000}, 
NDSolveValue[{ode[N1], T0[e] == 0, T0[1 - e] == 0}, {T0}, {x, e,1 - e}, 
WorkingPrecision -> 35][[1]]] 
Plot[T[x], {x, 0, 1}]

Using Ulrich's code above, adjust MaxStepSize, MaxSteps and WorkingPrecision as needed. This is the code for N1=30000 on a 3.5 GHz machine. Note I set up the problem for arbitrary precision arithmetic (use rational numbers for constants) and then set working precision to 25:

th = Pi/4;
s0 = 1/100 + 7/1000 Sin[th]
h = 3/10

ode[N1_] := (x (1 - x))/(2 N1) D[T0[x], x, x] + 
   s0 x (1 - x) (x + h (1 - 2 x)) D[T0[x], x] == -Integrate[
     E^(-N1 s0 a (a + 2 h (1 - a))), {a, 0, x}]/
   Integrate[E^(-N1 s0 a (a + 2 h (1 - a))), {a, 0, 1}]

AbsoluteTiming[
T=With[{e=10^-10,N1=30000},NDSolveValue[{ode[N1],T0[e]==0,T0[1-e]==0},{T0},{x,e,1-e},
MaxStepSize->1/20000,
MaxSteps->1000000,
WorkingPrecision->25][[1]]];
]

Out[61]= {26.3717,Null}