function [f,F] = func(x,data)
 

f=0.;

t=0;
dt = 0.01;
Dper=0.00154;
dx=0.02;
d=Dper/dx^2;
t2=0;
tfdi=0;
tfad=0;

numcel=20;
numpoints = 50000;
numiter = 13*1000*100;

period=[3 4 5 10];



% Known parameters
    uinf = 1;
    winf = 1;
    rinf=1;
    sinf = 1;
    vr=0.1;
    vc = 0.35;        % v_c
    vs = 0.6;        % v_s
    vo = 0;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% Unknown parameters
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

% Gate u
tumas  = x(1);
tumen  = x(2);

% Gate w
twmas  = x(3);
twmen  = x(4);

% Gate r
trmas  = x(5);
trmen  = x(6);

% Gate s
tsmas  = x(7);
tsmen  = x(8);

% Current xfi
gfi  = x(9);
vv   = x(10);

% Current xso
gso  = x(11);

% Current xsi
gsi  = x(12);
ks1  = x(13);
ks2  = x(14);
vs1  = x(15);
vs2  = x(16);

% Current xto
gto  = x(17);


% Initial values
for kcel = 1:numcel
  v=zeros(kcel,1);
  u=ones(kcel,1);
  w=ones(kcel,1);
  r=zeros(kcel,1);
  s=ones(kcel,1);
  lap=zeros(kcel,1);
end

  cv=zeros(4,1);
  apd=zeros(4,1);

for j = 1:4;
    
   voltage = zeros(numiter,1);
   s_save = zeros(numiter,1);
   r_save = zeros(numiter,1);
   
tfad=0;
tfdi=0;
   
for i = 1:numiter
for kcel = 1:numcel  
    
  if (kcel==1)  
     lap(kcel)=v(2)-v(1);
  elseif (kcel==numcel)  
      lap(kcel)=-v(numcel)+v(numcel-1);
  else     
      lap(kcel)=v(kcel+1)-2*v(kcel)+v(kcel-1);
  end
  
end

for kcel = 1:numcel 
  if(kcel<=3)  
   if (mod(i,period(j)*10000) == 9990) v(kcel) = 0.55;end
  end

    % Step functions
    p = (v(kcel) > vc);
    q = (v(kcel) > vs);
        
    % Gate evolution
    du = (1-p)*(uinf-u(kcel))/tumas - p*u(kcel)/tumen;
    dw = (1-p)*(winf-w(kcel))/twmas - p*w(kcel)/twmen;
    dr = q*(rinf-r(kcel))./trmas-(1-q)*r(kcel)./trmen;
    ds = (1-p)*(sinf-s(kcel))/tsmas - p*s(kcel)/tsmen;
    u(kcel) = u(kcel) + dt*du;
    w(kcel) = w(kcel) + dt*dw;
    r(kcel) = r(kcel) + dt*dr;
    s(kcel) = s(kcel) + dt*ds;
    
    % Currents
    xfi = -gfi*u(kcel)*p*(v(kcel)-vc)*(vv-v(kcel));
    xso = gso*((v(kcel)-vo)/vc*(1-p)+p);
    xsi = -gsi*w(kcel)*p*(1+tanh((v(kcel)-vs1)*ks1))*(1+tanh((vs2-v(kcel))*ks2));
    xto = gto*r(kcel)*v(kcel)*s(kcel);
   
    vold=v(kcel);

   v(kcel) = vold - dt*(xfi + xso + xsi+ xto)+dt*d*lap(kcel);
        
    


  if(kcel==numcel/2)

   voltage(i) = v(kcel);
   s_save(i) = s(kcel);
   r_save(i) = r(kcel); 

     if (v(kcel)>vr && vold<vr)
     tfdi=i*dt-(v(kcel)-vr)*dt/(v(kcel)-vold); 
     end
     if (v(kcel)<vr && vold>vr)
     tfad=i*dt-((v(kcel)-vr)*dt)/(v(kcel)-vold); 
     apd(j)=tfad-tfdi;
     end
   end

   if(kcel==numcel/2+1)
     if (v(kcel)>vr && vold<vr) 
     t2=i*dt-(v(kcel)-vr)*dt/(v(kcel)-vold);  
     end
   end
   
end

t = t + dt;


end

  if (t2~=0 && t2~=tfdi && tfdi~=0)
   cv(j)= dx/(t2-tfdi);
  else
   t2;
   tfdi;
  end
 


volt_sel = voltage(12*100000+1:10:12*100000+numpoints);
r_sel =r_save(12*100000+1:10:12*100000+numpoints); 
s_sel =s_save(12*100000+1:10:12*100000+numpoints); 

for k=1:numpoints/10
     if j~=1
            if (k < 1500) && (k > 1000)
                 a=20 ;
            else
                 a=1 ;
            end
    else
     a=1;
     end
  f=f+a*(volt_sel(k)-data(k,j))^2;
end

figure(1);subplot(1,4,j);plot(volt_sel,'k');hold on;plot(data(:,j),'r');hold off; drawnow;

end


f=f+10*((apd(4)-260.)^2+(apd(3)-257.)^2+(apd(2)-250.)^2+(apd(1)-226.)^2) ;


f=f+10000000*((cv(4)-0.072)^2+(cv(3)-0.071)^2+(cv(2)-0.068)^2+10*(cv(1)-0.0523)^2) ;

apd
cv
f
x

