I am having real trouble in creating a fourth order Runge-Kutta method for a system of two simultaneous first order ordinary differential equations.

It is supposed to model the wake off an aircraft wing, in terms of it's x and y coordinates.

The wing has been split into 10 parts, i.e. 5 on each side of the wing. And the program should give 10 values of x and y for each time step (i.e. 10 lots of 500 time steps in this case)

Any help with this problem would be gratefully appreciated..My attempted matlab program is below, but i fear ive got lost....

clc

m=10;

m2=m/2;

dx=1/m;

NMAX=500;

nfreq=1;

t=0;

x(1)=-1+dx;

ndt=0;

h=1;

dt=1/25/m;

for J=2:m;

x(J)=x(J-1)+2*dx;

end

for J=1:m2;

y(J)=0;

g(J)=((1-((x(J)+dx)^2))^0.5)-((1-((x(J)-dx)^2))^0.5);

end

m2p1=m2+1;

for J=m2p1:m;

y(J)=0;

g(J)=-g(m+1-J);

end

U=0;

V=0;

XI=-0.9;

YI=0;

NDT=0;

T=0;

while NMAX>NDT

NDT=NDT+1;

T=T+dt;

for J=1:10;

if J~=I;

U1=U+g(J)*(YI-Y(J))/((XI-X(J))^2 + (YI-Y(J))^2);

V1=V-g(J)*(XI-X(J))/((XI-X(J))^2 + (YI-Y(J))^2);

D1X= dt*U1;

D1Y= dt*V1;

X2=XI+D1X/2;

Y2=YI+D1Y/2;

U2=U+g(J)*(Y2-Y(J))/((X2-X(J))^2 + (Y2-Y(J))^2);

V2=V-g(J)*(X2-X(J))/((X2-X(J))^2 + (Y2-Y(J))^2);

D2X= dt*U2;

D2Y= dt*V2;

X3=XI+D2X/2;

Y3=YI+D2Y/2;

U3=U+g(J)*(Y3-Y2)/((X3-X2)^2 + (Y3-Y2)^2);

V3=V-g(J)*(X3-X2)/((X3-X2)^2 + (Y3-Y2)^2);

D3X= dt*U3;

D3Y= dt*V3;

X4=XI+D3X;

Y4=YI+D3Y;

U4=U+g(J)*(Y4-Y3)/((X4-X3)^2 + (Y4-Y3)^2);

V4=V-g(J)*(X4-X3)/((X4-X3)^2 + (Y4-Y3)^2);

D4X= dt*U4;

D4Y= dt*V4;

X(J+1) = X(J)+(D1X+2*D2X+2*D3X+D4X)/6.0;

Y(J+1) = Y(J)+(D1Y+2*D2Y+2*D3Y+D4Y)/6.0;

X(I)=X(J+1);

Y(I)=Y(J+1);

end

end

for I=m2p1:m;

X(I)=-X(m+1-I);

Y(I)=-Y(m+1-I);

end

X(1:10)

Y(1:10)

end

g(1:10)