Second Order ODE with Runge Kutta

Hello I need to solve the Water Droplet Trajectory with 4th Runge Kutta Method. My code runs but not gives me a logical values. Help me please.

Here is my code ( Just look after line 126)

% DAMLACIK YÖRÜNGE HESABI %

clc,clear all

V_inf =3; % Serbest Akım Hızı a = 1 ; %Yarıçap

[x,y]=meshgrid([-4:0.01:4]);

%%%% Silindir içerisinde akım fonksiyonuna sıfır değeri atanması

for i=1:length(x); for k=1:length(x); if (x(i,k).^2+y(i,k).^2)^(1/2)<a; x(i,k)=0; y(i,k)=0; end end end

%%%% DÖNGÜ SONU %%%%

r=sqrt(x.^2+y.^2); % SİLİNDİR YARI ÇAPI theta=atan2(y,x); % TETA AÇISI

P=V_inf.*sin(theta).*r.*(1-(a^2./(r.^2))); % Akım fonksiyonu

n=4000; %%%%%% r=ones(1,n+1).*a; z=[0:2*pi/n:2*pi];

%stream line Plot

figure(1) contour(x,y,P,50) hold on polar(z,r,'-b') axis square title('Silindir Etrafında Akım Çizgileri') hold on

x=[-a*2:0.04:a*2]; [x]=meshgrid(x); y=x'; for i=1:length(x); for k=1:length(x); if sqrt(x(i,k).^2+y(i,k).^2)<a; x(i,k)=0; y(i,k)=0; end end end

r=sqrt(x.^2+y.^2); theta=atan2(y,x); V_xr=V_inf*cos(theta).*(1-a^2./(r.^2)); V_yr=-V_inf*sin(theta).*(1+a^2./(r.^2));

V_x=V_xr.*cos(theta)-V_yr.*sin(theta); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% V_y=V_yr.*sin(theta)+V_x.*cos(theta); for i=1:length(x); for k=1:length(x); if sqrt(x(i,k).^2+y(i,k).^2)<a; V_x(i,k)=0; V_y(i,k)=0; end end end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

v=1.5*10^-5; %kinematik viskozite R=0.00001; %Yarıçap Ro= 1.225; %kg/m^3 yoğunluk Ro_w= 1000; %kg/m^3 SIVININ YOGUNLUGU g=9.81; % %Yer çekimi ivmesi Alfa=10; % Hücum açısı

d=2*R;

%%%% Sürtünme Katsayısının REynold Sayısına Göre Hesaplanması

% 10*10^-6 <d <500*10^-6 için

Re=V_inf*d/v; %Reynold Sayısı,

if 0<Re<0.05 C_d=24/Re; %Sürtünme Katsayısı

elseif 0.05<Re<3.04726 C_d= (24.167 + 3.254*Re - 0.23564*R^2)/Re ;

elseif 3.04726<Re<82.230 C_d= (-28.339 + 38.969*Re +0.73204*Re^2 - 0.0005608*Re^3)/Re^2;

elseif 82.230<Re

    C_d= exp(ln(Re)*(5.3220062-1.937238*lnRe + 0.4372856*(lnRe)^2 + 0.001926421*ln(Re)^4));
end

%%%% YÖRÜNGE DENKLEMİNİN HESAPLANMASI %%%%

K= 3*C_d*Ro/(4*d*Ro_w);

%RUNGE KUTTA ÇÖZÜMÜ %f1 = dVpx/dt %f2 = dVpy/dt %f3 = dx/dt %f4 = dy/dt

f1= @(Vpx,Vpy,t) K.*( (V_x(1,1)-Vpx).^2+(V_y(1,1)-Vpy).^2).^1/2.*(V_x(1,1)-Vpx) + g.*sin(Alfa); %dVpx/dt f3= @(Vpx,Vpy,t) Vpx; % dx/dt

f2= @(Vpx,Vpy,t) K.*( (V_x(1,1)-Vpx).^2+(V_y(1,1)-Vpy).^2).^1/2.*(V_y(1,1)-Vpy) + g.*cos(Alfa); % d Vpy/dt f4= @(Vpx,Vpy,t) Vpy; %dy/dt

Vpx(1)=V_x(1,1); Vpy(1)=V_y(1,1); t(1)=0; h=0.01; xxx(i)=5; yyy(i)=1; a=1; for i=1:100

    % k'lar Vpx hızı için
    % m'ler Vpy hızı için
    % l'ler x konumu için
    % n'ler y konumu için kullanılacak.

t(i+1)=t(i)+h;

k1=f1(Vpx(i),Vpy(i),t(i));  
m1=f2(Vpx(i),Vpy(i),t(i)); 
l1=f3(Vpx(i),Vpy(i),t(i));
n1=f4(Vpx(i),Vpy(i),t(i));

k2=f1(Vpx(i)+h*k1/2,Vpy(i)+h*m1/2,t(i)+h/2);
m2=f2(Vpx(i)+h*k1/2,Vpy(i)+h*m1/2,t(i)+h/2);
l2=f3(Vpx(i)+h*k1/2,Vpy(i)+h*m1/2,t(i)+h/2);
n2=f4(Vpx(i)+h*k1/2,Vpy(i)+h*m1/2,t(i)+h/2);

k3=f1(Vpx(i)+h*k2/2,Vpy(i)+h*m2/2,t(i)+h/2);
m3=f2(Vpx(i)+h*k1/2,Vpy(i)+h*m1/2,t(i)+h/2); 
l3=f3(Vpx(i)+h*k1/2,Vpy(i)+h*m1/2,t(i)+h/2);
n3=f4(Vpx(i)+h*k1/2,Vpy(i)+h*m1/2,t(i)+h/2);

k4=f1(Vpx(i)+h*k3,Vpy(i)+h*m3,t(i)+h);
m4=f2(Vpx(i)+h*k3,Vpy(i)+h*m3,t(i)+h);  
l4=f3(Vpx(i)+h*k3,Vpy(i)+h*m3,t(i)+h);  
n4=f2(Vpx(i)+h*k3,Vpy(i)+h*m3,t(i)+h); 

Vpx(i+1)=Vpx(i)+ h/6*(k1+2*(k2+k3)+k4);
Vpy(i+1)=Vpy(i)+ h/6*(m1+2*(m2+m3)+m4);
xxx(i+1) = x(i)  + h/6*(l1+2*(l2+l3)+l4); 
yyy(i+1) = y(i)  + h/6*(n1+2*(n2+n3)+n4);

end

figure(2) plot(t,Vpx)

figure(3) plot(t,Vpy)

figure(4) plot(t,xxx)

figure(5) plot(t,yyy)

figure(6) plot(xxx,yyy)