How to numerically solve differential equations with singularities?

Question

I wish to solve a differential equation numerically. As an example, I tried solving Hypergeometric differential equation

Block[{a = 1.23, b = 2.87, c = 2.76}, 
 s = NDSolve[{z (1 - z) w''[z] + (c - (a + b + 1) z) w'[z]- a b w[z] == 0, w[0] == 1, w'[0] == (a b)/c}, w, {z, -10, 10}, 
   SolveDelayed -> True]]

With the above when if I try to get the solution I get the following

In[]: Block[{a = 1.23, b = 2.87, c = 2.76, z = -1.23}, {Hypergeometric2F1[a, b, c, z], w[z] /. s}]
Out[]: {0.360842, {0.360842}}

which is indeed correct and gives me the correct result for any value $z<1$. But when I try to evaluate the same for some point $z>1$ then it gives a warning that "the point lies outside the range of data in the interpolation".

From what I understand this may be due to the singular point $z=1$ of the differential equation that might be the problem, but I am not able to find a way to surpass this problem. What can I do to resolve this issue?

As an added question given a general differential equation with singular points is it possible to solve it so as to get numerical results? Taking the example of $_2F_1$, if I have to find its value for points $|z| >1$ then I must use the analytic continuation formula. Is it possible to get the result by solving the differential equation numerically?

Edit 1: I tried a bit more but still it doesn't match the one evaluated using Mathematica's internal function. The result is complex but solving differential equations gives a completely real result which is plainly wrong. There is a branch cut for $[1,\infty)$.

Answer 1

The term hypergeometric function points to its origin, the geometric series

(1-x)^-1 = Sum[x^n,{n,0,oo}

This series has a pole at x=1 and converges in the complex open unit circle.

Its no problem, to use this series in greater circles centered at points on the real axis

(1-(x-a)^-1 = (1+a-x )^-1 = (1+a)^-1 * (1- x/(1+a))^-1

except for a=1. This is simply due to the fact that the function 1/x has no power series representation. It is its own Laurent series with one term.

It automatically follows, that all power series share the fact that their coefficients approach 1 for n->oo, producing a single pole at z=1. one has has to avoid the pole at z=1 in any evaluation procedure.

The defining series of the Gaussian hypergeometric functions all have quotients of products that grow as n! for large n and so the coeffienctients converge to 1 producing at least the pole at 1.

But this is not all: The coefficient series will contain subseries of the series converging as 1/n^k. They are producing logarithms and their integrals, all having a branch cut on (1,oo) (conventionally, or somewhere else on a curve joining the points 1 -- oo in the complex plane)

Of course there are exceptions: if the first parameter is a negative integer, the series terminates and yields a polynomial. Other cases avoiding logarithms exist.

Example

(Coefficient[ Normal@Series[Hypergeometric2F1[2, 3, 4, x], 
{x, 1, 4}] , {(x - 1)^-1, Log[1 - x]}]) // FullSimplify

{-3, 6 (35 + x (-105 + x (126 + 5 x (-14 + 3 x))))}

Hypergeometric2F1[-3, 3, 4, x]

1/20 (20 - 45 x + 36 x^2 - 10 x^3)

Rule of thumb for non-analysts: If in the definition of a function in the tables a branch cut is indicated, its always a good idea to assume that the complex argument makes a step of a ratonal multiople q*pi over the cut and and the complex values differ by exp(+-i q p). The standard case of a logarithm produces a change of sign, in rational powers its a factor on the unit circle.

The orignal question: How to circumvent the pole at z1 in the complex plane probably has two answers: Integrate on a little half circle around z=1 in the upper or lower plane. Depending on explicit powers and logs at z==1 the reuslts differ by the count of phases aquired.

Lets see:

 f = With[{a = 1.23, b = 2.87, c = 2.76},
 NDSolveValue[{z (1 - z) w''[z] + (c - (a + b + 1) z) w'[z] -  a b w[z] == 0,  
        w[2] == 1, w'[2] == (a b)/c},
        w, {z, 1.2, 10} ]];

g = With[{a = 1.23, b = 2.87, c = 2.76}, 
 DSolveValue[{z (1 - z) w''[z] + (c - (a + b + 1) z) w'[z] -  a b w[z] == 0, 
       w[2] == 1, w'[2] == (a b)/c}, 
       w, {z, 1.2, 10} ]]


Function[{z}, ((8.46197 + 1.35955 I) Hypergeometric2F1[-0.53, 
    1.11, -0.76, z])/z^( 44/25) - (10.5438 - 1.56523 I) Hypergeometric2F1[1.23, 2.87, 2.76,  z]]


Plot[{f[x], g[x]}, {x, 1.2, 10},  PlotStyle -> {{Red, Thickness[0.02]}, {Black}}]

ComplexPlot3D[g[z], {z, -1 - I, 2 + I}]

Its obvious that there is a branch cut along the positive real line x>1 and without exact description such definitions are of value for the complex analysts and theoretical physicists only.

The reason: Such integrals along cuts arise typically as solutions of wave equations of the second order in time and space by Fourier transforms. The two different solutions on the cut correspond to the split of the space of solutions into time forward and time backward. These two classes have obviously not very much in common, algebrically as phyically. The forward mode is used in preparing a signal sent to the receiving community in the future, the set of backward modes is used to analyze arriving signals.