I think there's a bug that leads NDSolve to conclude that the system is not separable. I think you should report it and see if WRI can verify it. It's a fair amount of work to track it down, and there is a lot of nearly unreadable stuff to wade through on the way.

To check separability, one can evaluate the P' derivatives at only the initial values of Q and, vice versa, the Q' derivatives at only the initial values of P. In both cases, if the result is a vector of numbers, then the system is separable. Unfortunately, NDSolve checks P' by plugging in the initial values of P for the variables Q. In many cases, this won't be a problem, since the number of inputs is the same and the results are thrown away anyway. But in this case, plugging the initial values of P into the variables Q causes divide-by-zero errors. Consequently the results for P' are not numeric and NDSolve says the system is not separable.

Here's a fix. NDSolve evaluates a result, but it will be up to the OP to determine whether it is reasonable result. The fix, as far as I can tell, only affects testing the system during initialization of the method; it should not affect the computation of the solution. (The fix consists of switching q0 and p0 in the two commented lines.)

cs[ndstate_, pnf_, ppos_, qnf_, qpos_] := (* checks separability *) Module[{dir, y0, p0, pf0, q0, qf0, sd}, dir = 1; sd = ndstate["SolutionData"[dir]]; y0 = ndstate["SolutionVector"[dir]]; q0 = y0[[qpos]]; p0 = y0[[ppos]]; Quiet[ NDSolve`SetSolutionDataComponent[sd, "X", p0]; (* Evaluate Q' at P0 *) qf0 = NDSolve`EvaluateWithSolutionData[qnf, sd]; NDSolve`SetSolutionDataComponent[sd, "X", q0]; (* Evaluate P' at Q0 *) pf0 = NDSolve`EvaluateWithSolutionData[pnf, sd];]; If[! VectorQ[qf0, NumberQ] || ! VectorQ[pf0, NumberQ], NDSolve`NDSolveMessage[ndstate, "sprksep", "SymplecticPartitionedRungeKutta"]; Throw[$Failed];]; ] Block[{NDSolve`SPRKDump`CheckSeparability = cs}, NDSolve[eqn, var, {t, 0, 10}, Method -> method, WorkingPrecision -> 10, MaxStepSize -> 0.001, MaxSteps -> ∞] ] 

I think there's a bug in the internal function NDSolve`SPRKDump`CheckSeparability that leads NDSolve to conclude that the system is not separable. I think you should report it and see if WRI can verify it (they would probably appreciate a link to this Q&A). It's a fair amount of work to track it down, and there is a lot of nearly unreadable stuff to wade through on the way.

To check separability, one can evaluate the P' derivatives at only the initial values of Q and, vice versa, the Q' derivatives at only the initial values of P. In both cases, if the result is a vector of numbers, then the system is separable. Unfortunately, NDSolve checks P' by plugging in the initial values of P for the variables Q. In many cases, this won't be a problem, since the number of inputs is the same and the results are thrown away anyway. But in this case, plugging the initial values of P into the variables Q causes divide-by-zero errors. Consequently the results for P' are not numeric and NDSolve says the system is not separable.

Here's a fix. NDSolve evaluates a result, but it will be up to the OP to determine whether it is reasonable result. The fix, as far as I can tell, only affects testing the system during initialization of the method; it should not affect the computation of the solution. (The fix consists of switching q0 and p0 in the two commented lines, basically a two-character fix.)

cs[ndstate_, pnf_, ppos_, qnf_, qpos_] := (* checks separability *) Module[{dir, y0, p0, pf0, q0, qf0, sd}, dir = 1; sd = ndstate["SolutionData"[dir]]; y0 = ndstate["SolutionVector"[dir]]; q0 = y0[[qpos]]; p0 = y0[[ppos]]; Quiet[ NDSolve`SetSolutionDataComponent[sd, "X", p0]; (* Evaluate Q' at P0 *) qf0 = NDSolve`EvaluateWithSolutionData[qnf, sd]; NDSolve`SetSolutionDataComponent[sd, "X", q0]; (* Evaluate P' at Q0 *) pf0 = NDSolve`EvaluateWithSolutionData[pnf, sd];]; If[! VectorQ[qf0, NumberQ] || ! VectorQ[pf0, NumberQ], NDSolve`NDSolveMessage[ndstate, "sprksep", "SymplecticPartitionedRungeKutta"]; Throw[$Failed];]; ] Block[{NDSolve`SPRKDump`CheckSeparability = cs}, NDSolve[eqn, var, {t, 0, 10}, Method -> method, WorkingPrecision -> 10, MaxStepSize -> 0.001, MaxSteps -> ∞] ] 

If I haven't got myself confused, I think it's an internal bug. I think you should report it and see if WRI can verify it. It's a fair amount of work to track it down, and there is a lot of nearly unreadable stuff to wade through on the way.

To check separability, one can evaluate the P' derivatives at only the initial values of Q and, vice versa, the Q' derivatives at only the initial values of P. In both cases, if the result is a vector of numbers, then the system is separable. Unfortunately, NDSolve checks P' by plugging in the initial values of P for the variables Q. In many cases, this won't be a problem, since the results are thrown away anyway. But in this case, plugging the initial values of P into the variables Q causes divide-by-zero errors. Consequently the results for P' are not numeric and NDSolve says the system is not separable.

Here's a fix. NDSolve evaluates a result, but it will be up to the OP to determine whether it is reasonable result. The fix, as far as I can tell, only affects testing the system during initialization of the method; it should not affect the computation of the solution.

cs[ndstate_, pnf_, ppos_, qnf_, qpos_] := (* checks separability *) Module[{dir, y0, p0, pf0, q0, qf0, sd}, dir = 1; sd = ndstate["SolutionData"[dir]]; y0 = ndstate["SolutionVector"[dir]]; q0 = y0[[qpos]]; p0 = y0[[ppos]]; Quiet[ NDSolve`SetSolutionDataComponent[sd, "X", p0]; qf0 = NDSolve`EvaluateWithSolutionData[qnf, sd]; NDSolve`SetSolutionDataComponent[sd, "X", q0]; pf0 = NDSolve`EvaluateWithSolutionData[pnf, sd];]; If[! VectorQ[qf0, NumberQ] || ! VectorQ[pf0, NumberQ], NDSolve`NDSolveMessage[ndstate, "sprksep", "SymplecticPartitionedRungeKutta"]; Throw[$Failed];]; ] Block[{NDSolve`SPRKDump`CheckSeparability = cs}, NDSolve[eqn, var, {t, 0, 10}, Method -> method, WorkingPrecision -> 10, MaxStepSize -> 0.001,  MaxSteps -> ∞] ] 

I think there's a bug that leads NDSolve to conclude that the system is not separable. I think you should report it and see if WRI can verify it. It's a fair amount of work to track it down, and there is a lot of nearly unreadable stuff to wade through on the way.

To check separability, one can evaluate the P' derivatives at only the initial values of Q and, vice versa, the Q' derivatives at only the initial values of P. In both cases, if the result is a vector of numbers, then the system is separable. Unfortunately, NDSolve checks P' by plugging in the initial values of P for the variables Q. In many cases, this won't be a problem, since the number of inputs is the same and the results are thrown away anyway. But in this case, plugging the initial values of P into the variables Q causes divide-by-zero errors. Consequently the results for P' are not numeric and NDSolve says the system is not separable.

Here's a fix. NDSolve evaluates a result, but it will be up to the OP to determine whether it is reasonable result. The fix, as far as I can tell, only affects testing the system during initialization of the method; it should not affect the computation of the solution. (The fix consists of switching q0 and p0 in the two commented lines.)

cs[ndstate_, pnf_, ppos_, qnf_, qpos_] := (* checks separability *) Module[{dir, y0, p0, pf0, q0, qf0, sd}, dir = 1; sd = ndstate["SolutionData"[dir]]; y0 = ndstate["SolutionVector"[dir]]; q0 = y0[[qpos]]; p0 = y0[[ppos]]; Quiet[ NDSolve`SetSolutionDataComponent[sd, "X", p0]; (* Evaluate Q' at P0 *) qf0 = NDSolve`EvaluateWithSolutionData[qnf, sd]; NDSolve`SetSolutionDataComponent[sd, "X", q0]; (* Evaluate P' at Q0 *) pf0 = NDSolve`EvaluateWithSolutionData[pnf, sd];]; If[! VectorQ[qf0, NumberQ] || ! VectorQ[pf0, NumberQ], NDSolve`NDSolveMessage[ndstate, "sprksep", "SymplecticPartitionedRungeKutta"]; Throw[$Failed];]; ] Block[{NDSolve`SPRKDump`CheckSeparability = cs}, NDSolve[eqn, var, {t, 0, 10}, Method -> method, WorkingPrecision -> 10, MaxStepSize -> 0.001, MaxSteps -> ∞] ]