function drds = ClassicalEqns(s,r)
global E info keep
% codes the semiclassical ODEs for multidimensional tunnelling; you shouldn't need to touch this function.
x = r(1); y = r(2);
dxds = r(3); dyds = r(4);

part = find(keep(:,1)<=x & keep(:,2)>=x & keep(:,3)<=y & keep(:,4)>=y);
part(2:end) =  [];
u = info{part,2}; ue = feval(u,x,y);
ux = info{part,3}; uxe = feval(ux,x,y);
uy = info{part,4}; uye = feval(uy,x,y);
V = 1./ue;
Vx = -uxe./ue.^2;
Vy = -uye./ue.^2;

drds(1) = r(3);
drds(2) = r(4);
drds(3) = (Vx - dxds.*(Vx.*dxds+Vy.*dyds))/2./(V-E);
drds(4) = (Vy - dyds.*(Vx.*dxds+Vy.*dyds))/2./(V-E);

drds = drds(:);
