% Simple Euler particle pushing program
tspan= [0. 2.*pi*1.]';
nn=100;

dt=(tspan(end)-tspan(1))/nn;
t=tspan(1):dt:tspan(end); %Calculate row vector of times

%Initial Conditions [ x  y vx vy ]'
y0= [0. 0. 1. 0.]'

%Parameters
B= 1;
qoverm=1.;

%Set solution as matrix of zeros
y=zeros(size(y0,1),size(t,2));

%Copy initial conditions y0  to first timestep of solution y
y(:,1)=y0(:);

%Loop to advance each timestep
for i = 2:nn+1
  dt=t(i)-t(i-1);
  y(1,i)=y(1,i-1) + dt*y(3,i-1);
  y(2,i)=y(2,i-1) + dt*y(4,i-1);
  y(3,i)=y(3,i-1) + dt*(qoverm*B*y(4,i-1));
  y(4,i)=y(4,i-1) + dt*(-qoverm*B*y(3,i-1));
end

%Calculate error due to numerical errors
terr=( (y(1,1)-y(1,end))^2. + (y(1,2)-y(2,end))^2.)^0.5

%Plot trajectory of particle (x,y)
plot(y(1,:),y(2,:))
xlabel('x')
ylabel('y')
title('Trajectory of Larmor Motion')