function F = Diffur (x, y);
% oscillational equations
global lorder lw
global pp_c2 pp_N2 pp_g pp_ro
global omegao2
%
c2x = ppval (x,pp_c2);
N2x = ppval (x,pp_N2);
gx = ppval (x,pp_g);
rox = ppval (x,pp_ro);
% c2x = interp1 (ro, c2o, x, 'cubic');
% N2x = interp1 (ro, N2o, x, 'cubic');
% gx = interp1 (ro, go, x, 'cubic');
% rox = interp1 (ro, rhoo, x, 'cubic');
rx = 1/x; % rx - reverse radius
c2r = 1/c2x; % inverse squared sound speed
lwo=lw/omegao2;
K = omegao2 - N2x; % + rox;        % for l=0 only, full solution
E = N2x/gx;
H = gx*c2r - 2*rx;
T = lwo*rx*rx - c2r;
MM=[E K; T H];
F = MM*y;
%F = [K*y(2) + E*y(1); H*y(2) + T*y(1)];