clear
clf
global qm E B betam 
q=1.e-9;m=1.e-9;beta=0.;qm=q/m;betam=beta/m;
 
E=[0.,0.1,0.];
B=[0.,0.,0.1];
y0=[0.,0.,0.,0.,0.,0];
options=odeset('AbsTol',1.e-9,'RelTol',1.e-9,'MaxStep',10.);
[t,y]=ode45(@lorentzf,[0,100],y0,options);
hold on

%axis([-1.5,1.5,-1.5,1.5,-1.5,1.5])
axis equal
plot3(y(:,1),y(:,2),y(:,3),'r');
quiver3(y(250,1),y(250,2),y(250,3),E(1),E(2),E(3),100)
quiver3(y(250,1),y(250,2),y(250,3),B(1),B(2),B(3),100)
grid on

xlabel('x');ylabel('y');

function dy = lorentzf(t,y)
global qm E B betam
dy=zeros(6,1);
dy(1)=y(4);
dy(2)=y(5);
dy(3)=y(6);
v=sqrt(power(y(4),2)+power(y(5),2)+power(y(6),2));
dy(4)=qm*(y(5)*B(3)-y(6)*B(2)+E(1))-betam*v*y(4);
dy(5)=qm*(y(6)*B(1)-y(4)*B(3)+E(2))-betam*v*y(5);
dy(6)=qm*(y(4)*B(2)-y(5)*B(1)+E(3))-betam*v*y(6);
end

