Python leapfrog trajectory code

I'm a physics student, but newbie in programming. Years ago I've learned how to do write a C code for a leapfrog integrator for a particle in gravitational field, but the memory is vague, and what I'm doing now is the writing code in Python for a leapfrog integrator for a particle in a certain magnetic field. Somebody told me Boris algorithm is better for this kind of simulation, but I decided first I would experiment with what I learned before, leapfrog integrator. But difference of syntax in C and Python was too great (at least to me) so I could not just translate the code from the C to Python, and I had to write a new one. So I'm not sure the algorithm is correct.

my code is like this,

# -*- coding: utf-8 -*-
"""
Created on Thu Feb 15 19:00:55 2018

@author: Heptacle
"""

import numpy as np
import matplotlib.pyplot as plt 

q=1.6e-19 # unit charge
m=1.67e-27 # proton mass
xs=1 # x_star
B0=1 # maximum magnetic field
b=B0/xs

initial_position = np.array((0.1, 0,0))         # Initial position vector of Particle
initial_velocity = np.array((0, 0,0.1))           # Initial velocity vector of Particle

num_steps = 4000
time_values = np.linspace(0, 1000, num_steps)
dt = time_values[1] - time_values[0]

positions = np.zeros((num_steps, 3))   
positions[0] = initial_position

velocities = np.zeros((num_steps, 3))    
velocities[0] = initial_velocity


def acc(x,v):
    if np.abs(x[0])<=1:
        B=(0,b*x[0],0)
    elif x[0]>=1:
        B=(0,B0,0)
    else:
        B=(0,-B0,0)
    a=q*np.cross(v,B)/m
    return a

vh=np.zeros((num_steps, 3))
vh[0]=velocities[0]+acc(positions[0],velocities[0])*dt/2

accs = np.zeros((num_steps, 3))
accs[0] = acc(positions[0],velocities[0])

for i in range(num_steps - 1):
    positions[i+1]=positions[i]+dt*vh[i]
    vh[i+1]=vh[i]+dt*acc(positions[i+1],vh[i])
    velocities[i+1]=vh[i+1]-dt/2*acc(positions[i+1],velocities[i])





####################
##### PLOTTING #####
####################
x_vals = positions[:,0]
y_vals = positions[:,1]
z_vals = positions[:,2]

plt.figure()
plt.plot(x_vals, z_vals, color = "blue", label = "Particle trajectory")
plt.legend(loc = "upper right")
plt.title("Orbit Plots")
plt.xlim((min(x_vals), max(x_vals)))
plt.ylim((min(z_vals), max(z_vals)))
plt.xlabel("x position ")
plt.ylabel("z position ")

plt.show()

But it didn't work properly. result plot

z value increases inexorably. It seems like some computational problem, but I cannot find exact problem. Can somebody help me?

I am not accustomed to the version of the leap-frog algorithm you are using. However, I tested your code and I think that the culprits are the q and m variables. The exceedingly small values they take are very likely to cause numeric issues. Indeed, my python interpreter even warns about this:

RuntimeWarning: overflow encountered in divide
  a=q*np.cross(v,B)/m

This sort of issues is very common in numerical simulations, and it can be handedly solved by using the so-called reduced units (see e.g. here). In your example, setting q = 1 and m = 1 seems to yield realistic results.

edit: I would like to add that using reduced units isn't as easy as setting all constants to 1, since physical units of measurements cannot be chosen independently. For instance, in molecular simulations (my field), it is customary to set to 1 the particle diameter, the energy scale and the mass. All other quantities are then expressed in these units. In your case, setting q = 1 and m = 1 would change the values of B0 and dt.