%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dydt = optAderiv(t,y, flag, yi, xi, p)

eta = 20;
M = length(xi);
diff = zeros(p, 1);

for i=1:M
  zetai = makezeta(xi(i), p);
  e = yi(i) - y'*zetai;
  diff = diff + zetai*e;
end

dydt = eta*diff;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define the basis functions
function [zeta] = makezeta(x, p)

sigma = 2;
cmax = 2*pi;
c = 1:p; c = (c(:) - 1)*cmax/(p-1);  % Gaussian centers

zeta = zeros(p,1);
for i=1:p
  zeta(i) = exp(-sigma^2*(x-c(i))^2);    
end
sz = sum(zeta);
zeta = zeta / sz;