function [dtheta_dt] = ODE_func(t, theta, gamma, Bmap1D, alpha)

c = 2.998e8; % speed of light (m/s)
m = 9.11e-31; % electron mass (kg)
q = -1.6e-19; % electron charge (C)

% params
x = theta(1);
y = theta(2);
z = theta(3);
vx = theta(4);
vy = theta(5);
vz = theta(6);

% get field at params
[ Bx, By ] = getfield( z, x, y, Bmap1D, alpha );

% changes in params
dxdt = vx;
dydt = vy;
dzdt = vz;
%  beta = sqrt(vx^2 + vy^2 + vz^2) / c;
%  gamma = 1./sqrt(1-beta^2);
%  display([z, gamma, vx/c, vy/c, Bx, By]);
%  display([ Bx, By]);
qbymgamma = q / m / gamma;
dvxdt = qbymgamma * - By * vz;
dvydt = qbymgamma * Bx * vz;
dvzdt = qbymgamma * (By * vx - Bx * vy);
%  dvzdt = 0.;

dtheta_dt = [dxdt, dydt, dzdt, dvxdt, dvydt, dvzdt]';

end