function [f] =  reiss(fi, y, Q, M, D)
% Reissner-Nordstrom null geodesic
% ODE for phi coordinate for the null geodesics
% f is d u / d phi
% for null geodesics with imaginary eccentricities
% Q, M and D are parameters
% fi and y are variables.
f = sqrt(-Q^2*y.^4 +2*M*y.^3 -y.^2 +1/D^2);