Error in using eval and ode

This is my code,

*% main function *

%F=farads
%V=mV
%A=amps
%F=Farads
C=1.0e-6;       %micro F/cm^2
gbar_Na=120; %(micro A/mV)/cm^2
gbar_K=36;   %(micro A/mV)/cm^2
gbar_L=0.3;   %(micro A/mV)/cm^2
E_Na=45;     %mV
E_K=-82;     %mV
E_L=-59;     %mv
%V=Y(1);
%m=Y(2);
%h=Y(3);
%n=Y(4);
Y0=[-69.89825095976047, 0.05357073962900, 0.59257698454360, 0.31924676125340];
%Y0(1)=Y0(1)+10;
t=0:0.01:60;
options=odeset('RelTol',1.e-8);

[T, Y]=ode45(@dydt_HH,t,Y0,options, C,gbar_Na,gbar_K,gbar_L,E_Na,E_K,E_L,...
             'Ioft');
%JR figure``(1)
figure(1);
plot(T,Y(:,1)) % time course of V(t)
xlabel('Time (ms)')
ylabel('V (mV)')
figure(2);
plot(T,(Y(:,1)+70)./100,T,Y(:,2),T,Y(:,3),T,Y(:,4)); % time courses of V (scaled),m,h,and n
legend('scaled V(t)','m(t)','h(t)','n(t)');
xlabel('Time'); ylabel('V and gating vars');
figure(3); clf;
hold on;
plot(Y(:,1),Y(:,2)); % V-m phase plane projection
Vpts=(-100.001:.5:40);
minfpts=m_and_tau_m(Vpts);
plot(Vpts,minfpts,'black');
xlabel('V (mV)'); ylabel('m');
figure(4); clf;
hold on;
plot(Y(:,1),Y(:,4)); % V-n phase plane projection
ninfpts=n_and_tau_n(Vpts);
plot(Vpts,ninfpts,'black');
xlabel('V (mV)'); ylabel('n');
figure(5);
plot(Y(:,4),Y(:,3)); % n-h phase plane projection
xlabel('n'); ylabel('h');

Subfunction 1:

function dY=dydt_HH(t,Y,...
                    C,gbar_Na,gbar_K,gbar_L,E_Na,E_K,E_L,Ifun)
V=Y(1);
m=Y(2);
h=Y(3);
n=Y(4);
[m_inf, tau_m]=m_and_tau_m(V);
[h_inf, tau_h]=h_and_tau_h(V);
[n_inf, tau_n]=n_and_tau_n(V);
I=eval([Ifun '(t)']);
dY=zeros(4,1);
dY(1)=(-1/C)*( gbar_Na*m^3*h*(V-E_Na) + gbar_K*n^4*(V-E_K)...
              + gbar_L*(V-E_L) - I );
dY(2)=(1/tau_m)*( m_inf-m );
dY(3)=(1/tau_h)*( h_inf-h );
dY(4)=(1/tau_n)*( n_inf-n );

Subfunction 2:

function [m_inf, tau_m]= m_and_tau_m(V)
alpha_m=( (V+45)/10 )./( 1-exp( -(V+45)/10 ) );
beta_m=4*exp( -(V+70)/18 );
m_inf=alpha_m./(alpha_m + beta_m);
tau_m=1./(alpha_m + beta_m);

Subfunction 3:

function [h_inf, tau_h]= h_and_tau_h(V)
alpha_h=0.07*exp( -(V+70)/20 );
beta_h=1./( 1 + exp( -(V+40)/10 ) );
h_inf=alpha_h./(alpha_h + beta_h);
tau_h=1./(alpha_h + beta_h);

Subfunction 4:

function [n_inf, tau_n]= n_and_tau_n(V)
alpha_n=0.1*( (V+60)/10 )./( 1-exp( -(V+60)/10 ) );
beta_n=0.125*exp( -(V+70)/80 );
n_inf=alpha_n./(alpha_n + beta_n);
tau_n=1./(alpha_n + beta_n);

I am getting error as,

Error using eval
Undefined function 'Ioft' for input arguments of type 'double'.

Error in dydt_HH (line 12)
I=eval([Ifun '(t)']);

Error in odearguments (line 88)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 114)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in hh (line 28)
[T, Y]=ode45(@dydt_HH,t,Y0,options, C,gbar_Na,gbar_K,gbar_L,E_Na,E_K,E_L,...

please can somebody help me in this regard.

Thanks in advance.

[Geoff Hayes]

Vetri - the error message is indicating a problem with the function Ioft which is your last input to the ode45

 [T, Y]=ode45(@dydt_HH,t,Y0,options, C,gbar_Na,gbar_K,gbar_L,E_Na,E_K,E_L,...
             'Ioft');

What is Ioft? Is it a function that you have defined somewhere?