function out = MyODEFunc
out{1} = @init;
out{2} = @fun_eval;
out{3} = []; %@jacobian;
out{4} = []; %@jacobianp;
out{5} = []; %@hessians;
out{6} = []; %@hessiansp;
out{7} = [];
out{8} = [];
out{9} = [];

% --------------------------------------------------------------------------
function dydt = fun_eval(t, x, C)

    
% PM state
V=x(1,1);
n=x(2,1);


% PM parameters
% units: V = mV; Kca = uM; Cm = uF/cm^2; g = uS/cm^2; phi = 1/s 
 Vk=-85;Vca=120;Kca=0.5;Cm=1;
 gk=20;gca=20;gkca=8;phi=12;
 V1=-3;V2=30;V3=-20;V4=30;

% sigma;eps; fi = unitless
sigma=0.185;eps=0.01;Vp=400;Kp=0.3;
% alpha = aMol cm^2/nC (approximately Acell/(2*Faraday) in those units)
alpha=.2;


% PM functions
minf=.5*(1+tanh((V-V1)/V2));
ninf=.5*(1+tanh((V-V3)/V4));
tau= 1/(phi*cosh((V-V3)/(2*V4)));
I_ca=gca*minf*(V-Vca);
I_k=gk*n*(V-Vk);
I_kca=gkca*(V-Vk)*C^4/(C^4+Kca^4);
jin=-alpha*I_ca;


% PM equations
xdot(1,1)= (-I_k - I_ca - I_kca)/Cm;
xdot(2,1)= (ninf-n)/tau;

dydt=xdot;

% --------------------------------------------------------------------------
function [tspan,y0,options] = init
handles = feval(MyODEFunc);
y0=[0,0];
options = odeset('Jacobian',handles(3));
tspan = [0 10];