You have to take the scale of your constants into account. The Lipschitz constant L of the ODE system can be taken as abs(Q/m*B)=1.6180874e+07. For the RK4 method to work sensibly you need L*h between 1e-4 and 1e-2, for the larger choice this is h=1e-9. The upper bound has then to be adapted to be conform with the number of integration steps, below I used 10000 steps to see a significant portion of the solution curve.
Next correct the use of indices, if yDot and yDotDot are on position 3 and 4 in the output dudt of uprime, then y and yDot have also to be found at these positions in the state vector.
R = 0.012;
a = 0.015;
b = 0;
w = 0.2;
v0 = 455*10^3;
function dudt = uprim(t, u)
m = 9.1091*10^-31;
Q = -1.6021*10^-19;
E = 20.0;
B = 0.92*10^-4;
x = u(1);
xDot = u(2);
y = u(3);
yDot = u(4);
r = sqrt( (x - a).^2 + (y - b).^2 );
xDotDot = (Q/m)*(E + B*yDot*(1 - w*(r/R)));
yDotDot = (Q/m)*( - B*xDot*(1 - w*(r/R)));
dudt = zeros(size(u));
dudt(1) = xDot;
dudt(2) = xDotDot;
dudt(3) = yDot;
dudt(4) = yDotDot;
end
h = 1e-9;
t = 0:h:1e-5;
u = zeros(6, length(t));
u(1, 1) = a - R;
u(2, 1) = v0;
u(3, 1) = 0;
u(4, 1) = 0;
u(5, 1) = 0;
u(6, 1) = 0;
for i = 1:(length(t) - 1)
k1 = uprim(t(i) , u(:, i) );
k2 = uprim(t(i) + 0.5*h, u(:, i) + 0.5*h*k1);
k3 = uprim(t(i) + 0.5*h, u(:, i) + 0.5*h*k2);
k4 = uprim(t(i) + h, u(:, i) + h*k3);
u(:, i + 1) = u(:, i) + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
end
figure(1);
plot(t, u(1,:), t, u(3,:));
figure(2);
plot(u(1,:),u(3,:));
The plot contains left the x and y components as functions of t and right the curve (x,y).
Let's call $\omega=\frac{qB_0}m$. Then your equations of motion are:
$$\begin{align}\frac{dv_x}{dt}&=0\\\frac{dv_y}{dt}&=\omega v_z\\\frac{dv_z}{dt}&=-\omega v_y\end{align}$$
Note the sign on the last equation. The initial conditions for velocity can be written as $v_x(0)=v_0$, $v_y(0)=v_0$, and $v_z(0)=0$. The first equation gives $v_x(t)=v_0$, so $x(t)=v_0t$.
Take the derivative of the second equation, and use $\frac{d v_z}{dt}$ from the third to get $$\frac{d^2v_y}{dt^2}+\omega^2 v_y=0$$
The solution is $$v_y(t)=A\sin(\omega t)+B\cos(\omega t)$$
To find out $A$ and $B$ you need the initial conditions. First
$$v_y(0)=B=v_0$$ Then: $$\frac{d v_y}{dt}(t)=\omega A\cos(\omega t)-v_0\omega\sin(\omega t)$$ At $t=0$ you have $$\frac{dv_y}{dt}(0)=\omega A=\omega v_z(0)=0$$
So you get $$v_y(t)=v_0\cos(\omega t)$$ and $$v_z(t)=-v_0\sin(\omega t)$$
Now integrate with respect to time, to get the positions. Using the initial conditions will shift the center for $z$
Best Answer
Seems to me you forgot the third equation involving the electric field:
$$ma_{x}=qv_{y}B \\ ma_{y}=-qv_{x}B \\ ma_z = -qE$$
Also, note that
$$\left(v_{x}^2+v_{y}^2+v_{z}^2\right)'=2(v_x a_x + v_y a_y + v_z a_z) = 0$$
In other words, extremal speeds are attained when velocity and acceleration are perpendicular. Thus,
$$v_x v_y B - v_x v_y B - v_z E = - v_z E = 0$$
Hence when the $z$-component of the velocity is $0$, you should be in an extremum of the speed.
However, your $z$-component of speed is not necessarily zero. In fact, what will happen is that the particle will be accelerated in the electric field vertically, while the magnetic field will induce a circular motion. Hence, the particle will move on a helicoidal trajectory. The minimum speed is therefore at the beginning of the motion, but this is not a stationary point of the equations.