Why can't I get a better fit?

I have been trying to fit a the data points to the function mentioned in the program using lsqcurvefit function by Levenberg Marquardt Algorithm. I tried numerous ways to get the fir but didn't find so. What can be the issue?

%Initializing the constants
clc; clear all; close all;
%For sodium current
gbar_Na=4;  %[mmho/cm^2]
gNaC=0.003; %[mmho/cm^2] 
ENa=50;     %[V}
%For slow inward currents
gs=0.09;    %[mmho/cm^2]
C=1;         %[microF/cm^2]
%Rate constants for gating variables a-in terms of alpha and b-in terms of
%beta
%For x1
C1_x1a=0.0005; C2_x1a=0.083; C3_x1a=50;  C4_x1a=0; C5_x1a=C4_x1a; C6_x1a=0.057; C7_x1a=1;
C1_x1b=0.0013; C2_x1b=-0.06; C3_x1b=20;  C4_x1b=0; C5_x1b=C4_x1b; C6_x1b=-0.04; C7_x1b=1;
%For m
C1_ma=0;       C2_ma=C1_ma;      C3_ma=47;   C4_ma=-1; C5_ma=47; C6_ma=-0.1;   C7_ma=C4_ma;
C1_mb=40;      C2_mb=-0.056; C3_mb=72;   C4_mb=0;  C5_mb=C4_mb;  C6_mb=C5_mb; C7_mb=C6_mb;
%For h
C1_ha=0.126;   C2_ha=-0.25;  C3_ha=77;   C4_ha=0; C5_ha=C4_ha; C6_ha=C5_ha;  C7_ha=C6_ha;
C1_hb=1.7;     C2_hb=0;      C3_hb=22.5; C4_hb=C2_hb;  C5_hb=C4_hb;  C6_hb=-0.082; C7_hb=1;
%For j
C1_ja=0.055;   C2_ja=-0.25;  C3_ja=78;   C4_ja=0; C5_ja=C4_ja;  C6_ja=-0.2;  C7_ja=1;
C1_jb=0.3;     C2_jb=0;      C3_jb=32;   C4_jb=0; C5_jb=C4_jb;  C6_jb=-0.1;  C7_jb=1;
%For d
C1_da=0.095;   C2_da=-0.01;  C3_da=-5;   C4_da=0; C5_da=0;  C6_da=-0.072;  C7_da=1;
C1_db=0.07;    C2_db=-0.017; C3_db=44;   C4_db=0; C5_db=0;  C6_db=0.05;    C7_db=1;
%For f
C1_fa=0.012;   C2_fa=-0.008; C3_fa=28;   C4_fa=0; C5_fa=0; C6_fa=0.15;    C7_fa=1;
C1_fb=0.0065;  C2_fb=-0.02;  C3_fb=30;   C4_fb=0; C5_fb=0; C6_fb=-0.2;    C7_fb=1;

%For stimulus
amp=-39.75;
n=1;
a=0;
c=400;
i=30;
T=1;

t1=[0:0.02:5000];
optn=odeset('MaxStep',T/3.5);
%Functions of alpha and beta for each gating parameter
%For x1           
alpha_x1=@(y)(C1_x1a*exp(C2_x1a*(y(1)+C3_x1a))+C4_x1a*(y(1)+C5_x1a))/(exp(C6_x1a*(y(1)+C3_x1a))+C7_x1a);
beta_x1=@(y)(C1_x1b*exp(C2_x1b*(y(1)+C3_x1b))+C4_x1b*(y(1)+C5_x1b))/(exp(C6_x1b*(y(1)+C3_x1b))+C7_x1b);
%For m
alpha_m=@(y)(C1_ma*exp(C2_ma*(y(1)+C3_ma))+C4_ma*(y(1)+C5_ma))/(exp(C6_ma*(y(1)+C3_ma))+C7_ma);
beta_m=@(y)(C1_mb*exp(C2_mb*(y(1)+C3_mb))+C4_mb*(y(1)+C5_mb))/(exp(C6_mb*(y(1)+C3_mb))+C7_mb);
%For h
alpha_h=@(y)(C1_ha*exp(C2_ha*(y(1)+C3_ha))+C4_ha*(y(1)+C5_ha))/(exp(C6_ha*(y(1)+C3_ha))+C7_ha);
beta_h=@(y)(C1_hb*exp(C2_hb*(y(1)+C3_hb))+C4_hb*(y(1)+C5_hb))/(exp(C6_hb*(y(1)+C3_hb))+C7_hb);
%For j
alpha_j=@(y)(C1_ja*exp(C2_ja*(y(1)+C3_ja))+C4_ja*(y(1)+C5_ja))/(exp(C6_ja*(y(1)+C3_ja))+C7_ja);
beta_j=@(y)(C1_jb*exp(C2_jb*(y(1)+C3_jb))+C4_jb*(y(1)+C5_jb))/(exp(C6_jb*(y(1)+C3_jb))+C7_jb);
%For d
alpha_d=@(y)(C1_da*exp(C2_da*(y(1)+C3_da))+C4_da*(y(1)+C5_da))/(exp(C6_da*(y(1)+C3_da))+C7_da);
beta_d=@(y)(C1_db*exp(C2_db*(y(1)+C3_db))+C4_db*(y(1)+C5_db))/(exp(C6_db*(y(1)+C3_db))+C7_db);
%For f
alpha_f=@(y)(C1_fa*exp(C2_fa*(y(1)+C3_fa))+C4_fa*(y(1)+C5_fa))/(exp(C6_fa*(y(1)+C3_fa))+C7_fa);
beta_f=@(y)(C1_fb*exp(C2_fb*(y(1)+C3_fb))+C4_fb*(y(1)+C5_fb))/(exp(C6_fb*(y(1)+C3_fb))+C7_fb);

%For Ik1 and Ix1bar
Ik1=@(y)0.35*(4*(exp(0.04*(y(1)+85))-1)/(exp(0.08*(y(1)+53))+exp(0.04*(y(1)+53)))+0.2*(y(1)+23)/(1-exp(-0.04*(y(1)+23))));
Ix1bar=@(y)0.8*(exp(0.04*(y(1)+77))-1)/exp(0.04*(y(1)+35));
%y(1)=Vm  y(2)=[Ca}l  y(3)=x1  y(4)=m  y(5)=h y(6)=j y(7)=d y(8)=f
%Defining the overall differential equation
F=@(t,y)[-(1/C)*(Ik1(y(1))+(Ix1bar(y(1))*y(3))+((gbar_Na*y(4)^3*y(5)*y(6)+gNaC)*(y(1)-ENa))+(gs*y(7)*y(8)*(y(1)+82.3+13.0287*log(y(2))))+stim1(t, amp, n, a, c , i, T));(-10^-7*(gs*y(7)*y(8)*(y(1)+82.3+13.0287*log(y(2))))+0.07*(10^-7-y(2)));(alpha_x1(y(1))*(1-y(3))-beta_x1(y(1))*y(3));(alpha_m(y(1))*(1-y(4))-beta_m(y(1))*y(4));(alpha_h(y(1))*(1-y(5))-beta_h(y(1))*y(5));(alpha_j(y(1))*(1-y(6))-beta_j(y(1))*y(6));(alpha_d(y(1))*(1-y(7))-beta_d(y(1))*y(7));(alpha_f(y(1))*(1-y(8))-beta_f(y(1))*y(8))];
[T1 Y]=ode15s(F,t1,[-84.62 1.8*10^-7 0.0074 0.0110 0.99 0.97 0.003 1],optn);

xdata=T1;
ydata=Y(:,1);

f=fit(xdata,ydata,'fourier1');
figure(1)
plot(xdata,ydata);
w0=f.w;
%figure(5)
%plot(f,xdata,ydata);
Ns=8;
N=length(ydata);
Z=[ones(N,1)];
for i=1:1:Ns
Z=[Z,cos(i*w0*t1'),sin(i*w0*t1')];
end
a=inv(Z'*Z)*(Z'*ydata);
Ak=[]; A0=a(1); Bk=[];
for i=2:2:length(a)-1
    Ak=[Ak;a(i)];
end
for i=3:2:length(a)
   Bk=[Bk;a(i)];
end
mag=[]; ph=[];
for i=1:1:length(Ak)
    m=sqrt((Ak(i)^2)+(Bk(i)^2));
    p=atan(Bk(i)/Ak(i));
    mag=[mag,m];
    ph=[ph,p];
end

f1=@(x,xdata)(x(1)+x(2)*cos(x(3)*xdata+x(4))+x(5)*cos(2*x(6)*xdata+x(7))+x(8)*cos(3*x(9)*xdata+x(10))+x(11)*cos(4*x(12)*xdata+x(13))+x(14)*cos(5*x(15)*xdata+x(16))+x(17)*cos(6*x(18)*xdata+x(19))+x(20)*cos(7*x(21)*xdata+x(22))+x(23)*cos(8*x(24)*xdata+x(25)));
x0=[A0 mag(1) w0 ph(1) mag(2) w0 ph(2) mag(3) w0 ph(3) mag(4) w0 ph(4) mag(5) w0 ph(5) mag(6) w0 ph(6) mag(7) w0 ph(7) mag(8) w0 ph(8)];



times=xdata(1):0.02:xdata(end);
options=optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
%options=optimoptions('lsqcurvefit','display','iter');
lb=[]; ub=[];

x=lsqcurvefit(f1,x0,xdata,ydata,lb,ub,options);


y1=f1(x,times);

figure(2)
plot(xdata,ydata,'ro',times,y1,'b');
xlabel('Time (ms)');
ylabel('Voltage');

The red line is the points plotted by getting the solution via using ode15s and the blue line is the solution from the function. KIndly tell me what can I do about it?