Please help me . ODE problem

I was trying to solve 4 first order differential equations, but there was something wrong. My equations :

dq3x/dt = z 
dq3y/dt =u
dz/dt=(-(q3x+5)/((q3x+5)^2+(q3y-10)^2)^(3/2)) -((q3x-5)/((q3x-5)^2+(q3y-10)^2)^(3/2))
du/dt=(-(q3y-10)/((q3x+5)^2+(q3y-10)^2)^(3/2)) -((q3y-5)/((q3x-5)^2+(q3y-10)^2)^(3/2))

And my code is given in the following

clc;                                               % Clears the screen
clear all;
%defined the position of 2 protons

h=0.01;
% interval time
t=0:h:10;
 
% initial conditions
q3x(1)=40;
q3y(1)=1;
z(1)=0;
u(1)=3;
%Function

F = @(t,q3x,q3y) z;
G = @(t,q3x,q3y) (-(q3x+5)/((q3x+5)^2+(q3y-10)^2)^(3/2)) -((q3x-5)/((q3x-5)^2+(q3y-10)^2)^(3/2));
P = @(t,q3x,q3y) u;
Q = @(t,q3x,q3y) (-(q3y-10)/((q3x+5)^2+(q3y-10)^2)^(3/2)) -((q3y-5)/((q3x-5)^2+(q3y-10)^2)^(3/2));

for i=1:(length(t)-1)
    k_1 = F(t(i),q3x(i),q3y(i));
    L_1 = G(t(i),q3x(i),q3y(i));
    M_1 = P(t(i),q3x(i),q3y(i));
    N_1 = Q(t(i),q3x(i),q3y(i));
    
    
    k_2 = F(t(i)+0.5*h,q3x(i)+0.5*h*k_1,q3y(i)+0.5*h*L_1);
    L_2 = G(t(i)+0.5*h,q3x(i)+0.5*h*k_1,q3y(i)+0.5*h*L_1);
    M_2 = P(t(i)+0.5*h,q3x(i)+0.5*h*M_1,q3y(i)+0.5*h*N_1);
    N_2 = Q(t(i)+0.5*h,q3x(i)+0.5*h*M_1,q3y(i)+0.5*h*N_1);
    
    k_3 = F((t(i)+0.5*h),(q3x(i)+0.5*h*k_2),(q3y(i)+0.5*h*L_2));
    L_3 = G((t(i)+0.5*h),(q3x(i)+0.5*h*k_2),(q3y(i)+0.5*h*L_2));
    M_3 = P((t(i)+0.5*h),(q3x(i)+0.5*h*M_2),(q3y(i)+0.5*h*N_2));
    N_3 = Q((t(i)+0.5*h),(q3x(i)+0.5*h*M_2),(q3y(i)+0.5*h*N_2));
    
    
    k_4 = F((t(i)+h),(q3x(i)+k_3*h),(q3y(i)+L_3*h));
    L_4 = G((t(i)+h),(q3x(i)+k_3*h),(q3y(i)+L_3*h));
    M_4 = P((t(i)+h),(q3x(i)+M_3*h),(q3y(i)+N_3*h));
    N_4 = Q((t(i)+h),(q3x(i)+M_3*h),(q3y(i)+N_3*h));
    
    q3x(i+1) = q3x(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  % main equation
    z(i+1) = z(i) + (1/6)*(L_1+2*L_2+2*L_3+L_4)*h;
    q3y(i+1) = q3y(i) + (1/6)*(M_1+2*M_2+2*M_3+M_4)*h;
    u(i+1)   = u(i) + (1/6)*(N_1+2*N_2+2*N_3+k_4)*h;
end    

[darova]

What is this for?

[darova]

Something is wrong with F and P functions. They always return 0 value

[Ameer Hamza]

Is it necessary to write your own numerical solver? Otherwise, you can use ode45().