ODE23s, Differential equation, problem solving

Hello everybody. I have a serious problem solving the diff. equation, expressed below. Any help;

function [kap, alf0, Pstar, Astar]= nofilter()

imu=sqrt(-1);

rho=1000;RR=0.3;

Jb=0.1; %diastato

CC=0.1;%diastato

chr=0.1;

cht=0.1;

U=0.5;

span=0.34/2;

alam1=90;

alam2=90;

chm=0.5*(chr+cht);

span3d=2*span;

area=2*(chm*span);

AR=(span3d^2)/area;

AR=10;

B=14.75; %adiastato

BB=B*rho*U*0.34*chr*chr*RR; %diastato

pic=0.25;

II=Jb/pi*rho*area*chr*chr*chr; %adiastato

C=CC/1/2*pi*rho*area*chr*U^2; %adiastato

a=span3d*RR^2/chr*area;

b=RR^2/chr^2;

e=RR/chr;

d=RR*span3d/area;

alf0=linspace(1, 60, 2);

kap=linspace(0.025,0.2,2);

for i=1:length(alf0)

for j=1:length(kap)

om(j)=pi*kap(j)*2*U/chr; % reduced velocity k, Huxham, k=pi*khuxham

T(j)=2*pi/om(j);

alf0r(i)=alf0(i)*pi/180;

theo(j)=besselh(1,2,kap(j))./(besselh(1,2,kap(j))+imu*besselh(0,2,kap(j))); %theo function, lift deficiency factor

Re=U*chr*10^6; %reynolds

CD0=0.075/(log10(Re)+2)^2; %skin friction

fprintf('Calculating for %d\n', i, j);

[t,y]=ode23s(@(t,y) nofilter(t,y, i, j), [0 , 100],[0, 0]);

[Astar(i,j), Pstar(i,j)] = processsignal(t,y);

Pstar(i,j) = 0.5*BB.*(Pstar(i,j))./0.5*rho* 2*RR.*(Astar(i,j)).*cos(((Astar(i,j)))).*U^3;

end

end

function dYdt=nofilter(t,Y, i, j)

thi(i,j)=alf0r(i).*cos(om(j)*t);

k1=8*II+2*a*cos(thi(i,j)).*(cos(Y(1))-Y(1)*sin(Y(1)));

dYdt=[Y(2);(1/k1)*((Y(2)).*(-2*B+theo(j).*(AR/AR+2).*b.*cos(thi(i,j)).*(-cos(Y(1)+Y(1).*sin(Y(1)))))+...

(Y(2).^2)*2*a*cos(thi(i,j))*(2*sin(Y(1))*Y(1)+Y(1)*cos(Y(1)))-Y(1)*C+...

theo(j)*(AR/AR+2)*e*(thi(i,j).*cos(thi(i,j))+(3/4-pic)*2*thi(i,j).*cos(thi(i,j)))+...

d*(thi(i,j)+(1/2-pic)*2*thi(i,j).*cos(thi(i,j)))-e*(1/pi)*sin(thi(i,j))*(CD0+1.6*(thi-atan(2*e*Y(2).*(cos(Y(1))-Y(1).*sin(Y(1))))).^2))];

end

end

function [Aout,Pstar] = processsignal(t,y)

% process signal

Y = y(:,1);

Ydot = y(:,2);

k = 2:length(Y)-1;

Kmax = find( Y(k)>Y(k-1) & Y(k)>Y(k+1) );

Kmax = Kmax+1; %correct index

Kmin = find( Y(k)<Y(k-1) & Y(k)<Y(k+1) );

Kmin = Kmin+1; %correct index

ti_max = 0.5*(t(Kmax(1:end-1))+t(Kmax(2:end)));

ti_min = 0.5*(t(Kmin(1:end-1))+t(Kmin(2:end)));

N = 5; % number of final cycles for averagingaw

Aout = 0.5*(mean(Y(Kmax(end-N,end))) - mean(Y(Kmin(end-N,end))));

subplot(2,1,1), plot(t,Y,t(Kmax),Y(Kmax),'ob',t(Kmin),Y(Kmin),'ob')

subplot(2,1,2), plot(t(Kmax),Y(Kmax),'o',t(Kmin),abs(Y(Kmin)),'o')

Jint = find(t>=t(Kmax(end-N)) & t<t(Kmax(end)));

Pstar = (1/4)*trapz(t(Jint),Ydot(Jint).^2); % integral,1/4

end