%%%%%%%% These are the default parameter values. Change for parameter estimation. 

k1 = 3/40; 
k2 = 1/25; 
k3 = 1/10; 
E1 = -280/3; 
Ena = 40; 
Edag = -80; 
Estar = -15; 
Fh = 1/2; 
Fn = 1/270; 
Gna = 100/3; 
g21 = -2; 
g22 = -9;

%%% Ensure the parameters are put in the "params" vector in this order (so CariqEq.m gets the correct parameter values)

params = [k1,k2,k3,E1,Ena,Edag,Estar,Fh,Fn,Gna,g21,g22];

%%%%%%%%

%%%%%%%% Do not change time range and initial conditions

tspan=[linspace(0,300,150) linspace(301,900,55)];

ICs = [-10 1 0];
%%%%%%%%


% Solve the model
solver_options = odeset('RelTol',1e-5); % Set relative tolerance to something small (1e-5)

[t,sol] = ode15s(@(t,sol) CaricEq(t,sol,params),tspan,ICs, solver_options);

% Output 
Eout = sol(:,1);
hout = sol(:,2);
nout = sol(:,3);


% Plot all the variables (E,h,n)

figure()
plot(t,sol(:,1),'r','LineWidth',2);

figure()
plot(t,sol(:,2:3),'b','LineWidth',2);


