[MATLAB数学相关] ode45求解Jiles-Atherton模型出现问题，求助！！

[zdx2001]

程序是参考另一个磁滞模型修改的，但应用JA模型始终图形不对，检查了几遍也不清楚问题出在哪，请各位高手，大侠帮忙一看，小弟谢过！！

代码如下：

function hysteresis
tspan=0:0.0015:10;
y0=0;
[t,y]=ode45(@odefun,tspan,y0);
A=500;
H=A*sin(2.3*t);
y;
H;
figure
plot(H,y)
xlabel('H')
ylabel('M')
title('Hysteresis curves')
size(y)
size(H)
function dydt=odefun(t,y)
alpha=8e-5;
c=0.55;
Ms=1.47e6;
k=60;
a=40;
A=500;
H=A*sin(2.3*t);
He=H+alpha*y;
if (He/a)<0.01
Man=Ms*((He/a)/3-((He/a)^3)/45+2*((He/a)^5)/945);
DMan=(1/a)*Ms*(1/3-((He/a)^2)/15+2*((He/a)^4)/189);

else
Man=Ms*(coth(He/a)-a/He);
DMan=Ms*((-1/a)*((csch(He/a))^2)+a/He^(-2));
end

Delt=sign(2.3*A*cos(2.3*t));
%DeltM=(1+Delt*sign(Man-y))/2;
dydt=2.3*A*cos(2.3*t)*((Man-y)+Delt*c*k*DMan)/(Delt*k-alpha*((Man-y)+Delt*c*k*DMan));
%dydt=2.3*A*sin(2.3*t)*(DeltM*(Man-y)+Delt*c*k*DMan)/(Delt*k-alpha*(DeltM*(Man-y)+Delt*c*k*DMan));