To elaborate on my comment, considerable progress can be made as follows. First, the constant of the motion given in my comment is derived by
system = Map[Factor, system, Infinity] /. B0 q/m -> b
(* {Derivative[2][x][t] == b (1 + k x[t]) Derivative[1][y][t],
Derivative[2][y][t] == -b (1 + k x[t]) Derivative[1][x][t]} *)
to combine multiple constants. Then,
Simplify[Subtract @@@ system.{x'[t], y'[t]}];
Simplify[2 Integrate[%, t]];
% == (% /. t -> 0 /. {Derivative[1][x][0] -> 10000, Derivative[1][y][0] -> 0})
(* Derivative[1][x][t]^2 + Derivative[1][y][t]^2 == 100000000 *)
Then, as noted by Alexei Boulbitch, the two ODEs in system are readily combined to give
system[[1]] /. y'[t] -> Sqrt[100000000 - Derivative[1][x][t]^2]
(* Derivative[2][x][t] == b (1 + k x[t]) Sqrt[100000000 - Derivative[1][x][t]^2] *)
However, at this point we shift the zero point of x[t] by 1/k and only then solve the resulting equation. (Of course, the initial condition on x[0] also must be shifted to 1/k.)
Flatten@DSolve[Derivative[2][x][t] == b k x[t] Sqrt[100000000 - Derivative[1][x][t]^2],
x[t], t]
(* x[t] -> InverseFunction[-((2 I EllipticF[I ArcSinh[Sqrt[(b k)/(20000 + 2 C[1])] #1],
(10000 + C[1])/(-10000 + C[1])] Sqrt[2 + (b k #1^2)/(-10000 + C[1])]
Sqrt[1 + (b k #1^2)/(20000 + 2 C[1])])/
(Sqrt[(b k)/(10000 + C[1])] Sqrt[400000000 - 4 C[1]^2 - 4 b k C[1] #1^2
- b^2 k^2 #1^4])) &][t + C[2]] *)
plus a second solution with is the negative of the first. They are much simpler expressions.
I should add that the two constants C can be evaluated after inverting the solution above, applying the boundary conditions, and performing a bit of algebra.
Alternative Derivation with Constants Determined
From my perspective, it is simpler and more instructive to solve the ODE using Integrate twice, each time evaluating the resulting constant of integration. Begin with
eq = (#/y'[t]) & /@ Simplify[system[[1]] /. x[t] -> x[t] - 1/k] /.
y'[t] -> Sqrt[100000000 - Derivative[1][x][t]^2]
(* b k x[t] == Derivative[2][x][t]/Sqrt[100000000 - Derivative[1][x][t]^2] *)
Note that the zero point of x[t] is shifted by 1/k, as above. Doing so enormously simplifies subsequent results. Next Integrate eq and add a constant of integration (equivalent to C[1] above).
i1 = Integrate[# x'[t], t] & /@ eq;
i1[[1]] = i1[[1]] + c1; i1
(* c1 + 1/2 b k x[t]^2 == -Sqrt[100000000 - Derivative[1][x][t]^2] *)
Now, apply the boundary conditions to determine and eliminate c1.
Flatten@Solve[i1, c1] /. t -> 0 /. {x'[0] -> 10000, x[0] -> 1/k}
(* {c1 -> -(b/(2 k))} *)
i1 = i1 /. %
(* -(b/(2 k)) + 1/2 b k x[t]^2 == -Sqrt[100000000 - Derivative[1][x][t]^2] *)
Next, Solve this equation for x'[t], so that Integrate can again be applied. (For now, consider only the first solution from Solve.)
eq2 = Equal @@ (Solve[i1, x'[t]][[1, 1]]);
eq2 = eq2[[1]]/eq2[[2]] == 1
(* -((2 k Derivative[1][x][t])/
Sqrt[-b^2 + 400000000 k^2 + 2 b^2 k^2 x[t]^2 - b^2 k^4 x[t]^4]) == 1 *)
i2 = Reverse[Integrate[# , t] & /@ eq2];
i2[[2]] = i2[[2]] + c2; i2
(* t == c2 + (2 I k EllipticF[I ArcSinh[Sqrt[-((b k^2)/(b - 20000 k))] x[t]],
(b - 20000 k)/(b + 20000 k)] Sqrt[1 - (b k^2 x[t]^2)/(b - 20000 k)]
Sqrt[1 - (b k^2 x[t]^2)/(b + 20000 k)])/(Sqrt[-((b k^2)/(b - 20000 k))]
Sqrt[-b^2 + 400000000 k^2 + 2 b^2 k^2 x[t]^2 - b^2 k^4 x[t]^4]) *)
A substantial simplification now can be obtained by using the identity,
Sqrt[1 - (b k^2 x[t]^2)/(b - 20000 k)] Sqrt[1 - (b k^2 x[t]^2)/(b + 20000 k)] ==
Sqrt@Expand[-(b - 20000 k - b k^2 x[t]^2) (b + 20000 k -
b k^2 x[t]^2)] Sqrt[-1/(b - 20000 k)] Sqrt[1/(b + 20000 k)]
Applying it yields
(* t == c2 + (2 I EllipticF[I ArcSinh[Sqrt[-(b/(b - 20000 k))] k x[t]],
(b - 20000 k)/(b + 20000 k)])/Sqrt[b (b + 20000 k)] *)
Finally, determine and eliminate c2, and shift the zero point of x[t] back.
Simplify[Flatten@Solve[i2 /. t -> 0 /. x[0] -> 1/k, c2], k > 0 && b > 0]
(* {c2 -> -((2 I EllipticF[I ArcSinh[Sqrt[-(b/(b - 20000 k))]], (b - 20000 k)/
(b + 20000 k)])/Sqrt[b (b + 20000 k)])} *)
i2 = Simplify[i2 /. % /. x[t] -> x[t] + 1/k]
(* t == -((2 I (EllipticF[I ArcSinh[Sqrt[-(b/(b - 20000 k))]],
(b - 20000 k)/(b + 20000 k)] -
EllipticF[I ArcSinh[Sqrt[-(b/(b - 20000 k))] (1 + k x[t])],
(b - 20000 k)/(b + 20000 k)]))/Sqrt[b (b + 20000 k)]) *)
(Note that the second solution from Solve, alluded to earlier, results in the negative of the expression just obtained.) These results, although cumbersome, do provide analytical expressions for the limits on x[t] and for its oscillation period.
xlim = Solve[b/(b - 20000 k) k^2 x[t]^2 == (b + 20000 k)/(b - 20000 k), x[t]]
(* {{x[t] -> -(Sqrt[b + 20000 k]/(Sqrt[b] k))}, {x[t] -> Sqrt[b + 20000 k]/(Sqrt[b] k)}} *)
period = Simplify[2 Simplify[Subtract @@ (i2[[2]] /. x[t] -> x[t] - 1/k /. %)],
k > 0 && b > 0]
(* -((8 I EllipticF[I ArcSinh[Sqrt[-((b + 20000 k)/(b - 20000 k))]],
(b - 20000 k)/(b + 20000 k)])/Sqrt[b (b + 20000 k)]) *)
To illustrate these results graphically, consider the case {b -> 1, k -> 1}.
b1k1 = i2 /. {b -> 1, k -> 1} /. x[t] -> x;
shift = % /. {b -> 1, k -> 1};
Show[ParametricPlot[{shift # + (2 Mod[#, 2] - 1) b1k1[[1]], x},
{x, -1 - Ceiling[#, 6] Sqrt[20001], -1 + Sqrt[20001]},
AspectRatio -> 1/GoldenRatio, AxesLabel -> {t, x}] & /@
Range[0, 5], PlotRange -> {{0, .2}, Automatic}]
Numerical Solution
The corresponding numerical solution can be obtained without difficulty.
system = {x''[t] == ((((B0 k) x[t] + B0) q)/m) y'[t],
y''[t] == -((((B0 k) x[t] + B0) q)/m) x'[t]} /. {B0 -> 1, q -> 1, m -> 1, k -> 1};
initialvalues2 = {x[0] == 0, y[0] == 0, x'[0] == 10^4, y'[0] == 0};
solution = Flatten@NDSolve[{system, initialvalues2}, {x[t], y[t], x'[t], y'[t]},
{t, 0, 2}];
Plotting x[t] /. solution gives a figure identical to that obtained from the analytical solution above. For completeness, all four numerical curves are given by
Plot[Evaluate[{x[t], y[t]} /. solution], {t, 0, .2}]
Plot[Evaluate[{x'[t], y'[t]} /. solution], {t, 0, .2}]
I would close by observing that, while an analytical solution indeed can be obtained, a numerical solution is easier to derive and probably is more useful.
Derivative[1][x][t]^2 + Derivative[1][y][t]^2is a constant of the motion. $\endgroup$