# How to specify initial condition including integral equations in the case of the dirac equations?

Hi I am trying to solve a system of differential equation with NDSolve. The problem is I have like $3$ (maybe $5$) boundary conditions, but only $2$ differential equations. The equations to solve are the Dirac equations in a hydrogen like potential. My code is:

    k := -1; n := 1; Z := 82; α := 1/137
u[r_] := (Z*α)/r; ζ := (1 + (Z*α)^2/(n - Abs[k] + Sqrt[k^2 - (Z*α)^2])^2)^(-1/2)
eqn := {g'[r] == -(k/r) g[r] + (ζ + 1 + u[r]) f[r],
f'[r] == k/r f[r] - (ζ - 1 + u[r]) g[r]}
bound := {g[10^(-20)] == 0, f[10^(-20)] == 0}

dirac = NDSolve[{eqn, bound}, {g, f}, {r, 10^(-10), 20},
Method -> {"ExplicitRungeKutta"}]

Plot[Evaluate[{g[r], f[r]} /. dirac], {r, 0, 15}, PlotRange -> {{0, 15}, {-0.2, 0.6}}]


If I run this code, the solutions will be the $0$ functions. But that is not what I want.

The other boundary condition is The normalization of those functions, but I don't know how to express it as a boundary condition:

    Integrate[Abs[g[r]]^2 + Abs[f[r]]^2, {r, 0, Infinity}] == 1


This condition somehow include the conditions:

    g[Infinity]==0, f[Infinity]==0


If I exchange one of my condition with the integral or the others, mathematica won't find a solution.

The solution of those equations should be 0 at r = 0 and they should vanish for r = ∞, but the area under the curve should not be 0.

I hope someone can help me soon.