function F = Diffur (x, y)
% oscillational equations
global c2o N2o go ro lorder
global omegao2
lw = lorder*(lorder+1);
c2x = interp1 (ro, c2o, x, 'cubic');
N2x = interp1 (ro, N2o, x, 'cubic');
gx = interp1 (ro, go, x, 'cubic');
rx = x;
K = omegao2 - N2x;
E = N2x./gx;
H = gx./c2x -2./rx;
T = (lw./((rx.^2)*omegao2)) - (1./c2x);

F = [K*y(2) + E*y(1); H*y(2) + T*y(1)];