Do not use the for loop here. Let vectorisation work for you.
Try this:
eta=0.25;
L=2.5;
a=0.1;
h=0.1;
b=6;
x=a:h:b;
A =@(x) (0<x).*1+(0<=x & x<=L).*(1+0.5*(1-cos(2*pi*x/L)).^1/2) + (L<x).*1;
A1=@(x) (gradient(A(x))./gradient(x));
f1=@(x,y) A1(x)+(1.33*pi*eta*y.^3)./(1+1.33*pi*eta*y.^3);
y(1)=01;
k1=h*f1(x,y);
k2=h*f1(x+0.5*h,y+.5*k1);
k3=h*f1(x+.5*h,y+.5*k2);
k4=h*f1(x+h,y+k3);
y = +(1/6)*(k1+2*k2+2*k3+k4);
plot(x,y)
.
Best Answer