[t,u,uprime,rpar,ipar]=fdeqc(t,u,uprime,rpar,ipar);

function [t,u,uprime,rpar,ipar]=fdeqc(t,u,uprime,rpar,ipar);
persistent r r3 rsq ; 

ipar_shape=size(ipar);ipar=reshape(ipar,1,[]);
rpar_shape=size(rpar);rpar=reshape(rpar,1,[]);
u_shape=size(u);u=reshape(u,1,[]);
uprime_shape=size(uprime);uprime=reshape(uprime,1,[]);
if isempty(r), r=0; end;
if isempty(rsq), rsq=0; end;
if isempty(r3), r3=0; end;
rsq = u(1).*u(1) + u(2).*u(2);
r = sqrt(rsq);
r3 = rsq.*r;
uprime(1) = u(3);
uprime(2) = u(4);
uprime(3) = -(u(1)./r3);
uprime(4) = -(u(2)./r3);
ipar_shape=zeros(ipar_shape);ipar_shape(:)=ipar(1:numel(ipar_shape));ipar=ipar_shape;
rpar_shape=zeros(rpar_shape);rpar_shape(:)=rpar(1:numel(rpar_shape));rpar=rpar_shape;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
uprime_shape=zeros(uprime_shape);uprime_shape(:)=uprime(1:numel(uprime_shape));uprime=uprime_shape;
end %subroutine fdeqc